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INTRODUCTION 


Objectives and Statement of the Problem 

Inducers are designed to increase the suction specific speed of 
centrifugal impellers in liquid rocket feed systems and are used to 
provide a small head rise sufficient to operate without cavitation. 
Typical inducer characteristics include high solidity (ratio of blade 
chord to blade spacing) , low aspect ratio (span squared to blade area) , 
and low flow coefficient (ratio of inlet axial velocity to blade tip 
speed) . The flow in these long and narrow passages is greatly 
influenced by the effects of turbulence and viscosity, resulting in 
large friction losses and introducing considerable three-dimensionality 
in the flow, thus making the prediction of the flow behavior extremely 
difficult. Secondary motions within the inducer are not confined to 
thin regions at the blade surface, but extend over the entire cross- 
section of the flow. 

The primary objective of the current study reported in this thesis 
has been to gather additional understanding of the complex inviscid and 
viscid effects on the inducer flow field, including the three components 
of mean velocity, turbulence intensities and turbulence stresses inside 
the passage. It is hoped that the knowledge gained from this investi- 
gation will serve the establishment of a theoretical model for the 
eventual analysis of the three-dimensional flow in inducers as well as 
other turbomachinery dominated by secondary fluid motions caused by 
viscosity and turbulence. Hence, the subjects addressed in this thesis 


are: 
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1) To develop measuring techniques using a triaxial hotwire probe 
to measure the three components of velocity, turbulence intensities and 
stresses within the rotating blade channel. 

2) Using these techniques, to carry out a complete flow survey of 
two axial stations within the blade passage. 

3) To perform a complete survey of the static pressure distribution 
on the blade suction and pressure surfaces. 

4) To perform a complete survey of limiting streamline angle on 
the blade suction and pressure surfaces. 

5) To predict the three-dimensional inviscid flow in the inducer 
by numerically solving the exact equations of motion and to study the 
possible methods of reducing the computation time required for the 
convergence to the solutions of these equations. 

6) To approximately predict the three-dimensional viscid flow in 
the inducer using dominant viscous terms in the exact solution equations. 

Previous Related Work 

The Department of Aerospace Engineering at The Pennsylvania State 
University has been conducting a systematic analytical and experimental 
investigation of flow behavior in axial flow inducers under NASA 
sponsorship since November 1963. A brief summary of previous 
theoretical and experimental results obtained from this investigation 
is given in this section, in addition to the applicable results of 
related research by other sources. 
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Analytical Investigations 

Because of the presence of large secondary flows caused by three- 
dimensional boundary layers and the complexity of the viscid equations 
of motion governing the inducer flow, very little theoretical analysis 
is available related to the prediction of the three-dimensional flow 
characteristics. Most of the design and analysis of the inducer fluid 
flow is based on conventional two-dimensional methods. 


Montgomery (Ref. 17) used the simplified radial equilibrium 


equation 


r 

3h 

8r = v 


CD 


in conjunction with arbitrary expressions for the losses to predict the 
exit head rise and flow coefficient for an 80° helical inducer. However, 


the use of arbitrary loss expressions does not provide any specific 
method of relating the loss distribution to a given inducer geometry or 
flow characteristic. 


An approximate solution using the simplified radial equilibrium 
equation in an integrated form has been obtained for a four -b laded 
inducer by Lakshminarayana (Ref. 9) . The basic assumptions in this 
analysis are the existence of fully developed turbulent flow, the shape 
of the radial and mainstream direction velocity profiles and the use of 
empirically determined friction loss coefficients. This analysis leads 
to a good prediction of the outlet absolute tangential velocity, but 
axial velocities are not predicted due to the lack of accurate information 
concerning the radial velocity profiles within the blade passage. 

The momentum integral equations valid for the inducer have been 
developed and programmed for numerical solution (Ref. 2) . The equations 
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take into account the changes in free stream velocity due to camber and 
incidence in addition to the interaction between the pressure and 
suction surface boundary layers, and is based on skin friction 
correlations for rotating boundary layers developed in Ref. 14. The 
numerical technique utilizes a fourth-order difference scheme. 

Preliminary results of the analysis are discussed in Ref, 2. 

An accurate knowledge of boundary layer characteristics and skin 
friction losses in a rotating channel is a prerequisite for the develop- 
ment of an acceptable theoretical model for the inducer flow. While a 
considerable amount of information on viscous flow in a non-rotating flow 
passage is available, there is no information available for the rotating 
case. A systematic study of the boundary layer on a rotating blade and 
inside a rotating channel has been undertaken by the Aerospace 
Engineering Department at The Pennsylvania State University. In the 
first phase of this program, the boundary layer on a simpler configura- 
tion (a single rotating helical blade of large chord length enclosed in 
an annulus) was studied. The results of this investigation are reported 
in Ref. 14. Consequently, this investigation has been extended to a 
four-bladed flat plate inducer in Ref. 2, where an attempt has been made 
to predict and measure the boundary layer characteristics inside the 
blade passage. Further study is to be made on the three-dimensional 
boundary layer characteristics within the inducer passage by utilizing 
the rotating hotwire anemometry techniques developed in this thesis. 
Information gained through this investigation should provide knowledge 
of the significant viscous effects within the blade passage and 
eventually lead to the exact prediction of the inducer flow. 
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A significant contribution toward the general solution of the 
equations governing the inducer flow is due to Cooper and Bosch (Ref. 4). 
This three-dimensional analysis employs an iterative numerical procedure 
to solve the inviscid equations of motion, expressed in finite -difference 
form, for a grid of points representing the channel between the blades. 
This method is discussed in greater detail in a later section of this 
thesis . 

Experimental Investigations 

Several experimenters have investigated inducers of varying 
geometry and inlet angle and have tested them in various fluids such 
as water, liquid hydrogen or nitrogen under a wide range of flow 
parameters (Refs. 16, 17, 18, 19, 21, 23, 24). In most cases, however, 
these studies deal only with cavitation performance, overall performance 
and efficiency. 

Acosta (Ref. 1) has studied the cavitating and non-cavitating 
performance of 78°, 81° and 84° helical flat plate inducers under 
various flow coefficients. His investigations at <j> = 0.070 under 
non-cavitating conditions observed a deterioration in the radial 
distribution of axial velocity and head rise at the inducer exit, with 
a backflow region near the hub and a sharp positive gradient in head rise 
near the tip. These results are similar to those obtained in Ref. 20 
for a three-bladed inducer at <J> = 0.065. Acosta attributes strong three- 
dimensional and viscid effects for the departure of the flow from design 
values based of the simplified radial equilibrium equation. He also 
noticed a decrease in inducer efficiency at large solidity, which can be 
attributed to the influence of blade blockage on flow characteristics 
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and an increase in viscous and turbulent mixing losses due to the 
decrease in channel width. These results are confirmed in Refs. 11 and 
20 . Soltis, Anderson and Sandercock (Ref. 23 ) were led to similar 
conclusions while investigating the non-cavitating performance of a 78° 
axial inducer under various flow coefficients. They derived the outlet 
axial velocity profile using experimental values of the total pressure 
and outlet flow angles in the simplified radial equilibrium equation. 
This analysis tends to establish that the flow is axisymmetric at small 
axial distances downstream of the trailing edge, since radial velocities 
are likely to be small and the wake diffusion in such inducers is very 
rapid. Similar observations are made by Mullan (Ref. 18), Meng and 
Moore (Ref. 16), Montgomery (Ref. 17), and Osborn (Ref. 19). 

The main conclusions of the various investigations described above 

are: 

1) The overall head rise coefficient increases, especially near 
the tip, when the operating flow coefficient decreases. 

2) The total head rise coefficient increases when the solidity of 
the blades is decreased. 

3) The radial distribution of outlet velocity tends to deteriorate 
when the flow coefficient is decreased. At low flow coefficients and 
for most inducer configurations, there is a large positive radial 
gradient in exit axial velocity with a backflow near the hub. 

Experimental investigations at The Pennsylvania State University 
have been conducted on a three-foot diameter inducer operated in air at 
a flow coefficient of (p = 0.065. The inducer has been operated with 
four, three, and two blades giving solidities at the tip of 2,86, 2.13, 
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and 1.43 respectively. The inducer test facility is shown in Fig. la, 
with design values of inducer inlet and outlet angles given in Fig. lb. 

A visualization study of the flow through the four-b laded inducer 
configuration is reported in Ref. 10. The flow near the blade surfaces, 
inside the rotating passage, downstream and upstream of the inducer is 
visualized by means of smoke, tufts, ammonia filament and lamp black 
techniques. The flow is found to be highly three-dimensional with 
appreciable radial velocity throughout the passage. Some of the major 
conclusions of this visualization study are: 

1) At or near design flow coefficient, no backflow is observed up- 
stream of the inducer. A separated region of the flow exists near the 
hub at the discharge of the inducer. 

2) The extent of the backflow increases considerably, both at 
inlet and at exit, for flow coefficients lower than the design value. 

3) The expected radial motions within the blade passage have been 
confirmed and appear to be quite strong at all the radii. 

4) The radial flows inside the blade boundary layer, when 
encountered by the annulus wall, tend to deflect toward the mid-passage 
and then radially inward. 

The qualitative nature of the velocity profiles, derived from 
visualization experiments (Ref. 10), indicate that the conventional 
practice of assuming the boundary layer is thin, two-dimensional and is 
a small perturbation of the inviscid flow is not valid in this case. 

The values of limiting streamline angles are found to be large. 

The flow measurement at several stations downstream of the blade 
row are reported in Refs. 9, 20 and 11 for four-, three- and two-bladed 
inducers respectively. The radial distribution of stagnation and static 
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pressure, axial and tangential velocity and flow angles of the absolute 
flow were measured at the several locations downstream using conventional 
and hotwire probes. Absolute flow measurements were also carried out 
at the exit of the three-bladed and two-bladed inducers and reported in 
Refs. 8 and 20 and Ref. 11 respectively. All inducers were tested at 
the same flow coefficient and Reynolds number. Major conclusions 
derived by comparison of these measurements with those of the four- 
fa laded inducer are: 

1) The performance of the inducer improves continuously with 
decrease in solidity, the two-bladed inducer showing substantial improve- 
ment over both the four- and three-bladed inducers. 

2) The static and stagnation head rise increases continuously, at 
all radii, with decrease in solidity. The radial gradient of stagnation 
head rise coefficient, 8^/Sr, is found to be almost constant from hub 
to tip for the two-bladed inducer, unlike that at the other solidities 
where a steep rise is observed neaT the tip. 

3) The downstream axial velocity profile is found to be similar, 
qualitatively, for all the inducers tested. The steep rise in axial 
velocity toward the tip observed in three- and four-bladed inducers is 
absent in the case of the two-bladed configuration, but the extent of 
separated zone (backflow) near the hub increases continuously with 
decrease in solidity. 

4) The radial distribution of tangential velocity shows a trend 
similar to the ^ distribution. However, the large values of 9V/3r 
observed in four- and three-bladed inducers are reduced substantially in 
the two-bladed configuration. 
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5) Hydraulic efficiency (local as well as overall) increases 
continuously with decrease in solidity. 

In order to understand the flow behavior inside the inducer passages, 
experimental investigation of the relative b lade-to-blade flow inside the 
inducer passages is undertaken using rotating pressure probes and the 
pressure transfer device. Measurements inside the four-b laded inducer 
are reported in Refs. 7 and 15 and for the three-bladed inducer in Refs. 

7 and 20. The measurement of relative flow near the trailing edge 
reveals the presence of a loss core located slightly inward from the 
tip. The mid-passage at this radius is found to have minimum relative 
stagnation pressure and hence maximum loss. The radial velocity inside 
the blade boundary layer, when encountered by the annulus wall, tends to 
deflect toward the mid-passage and then radially inward. These inter- 
action effects are responsible for the large losses observed experi- 
mentally. These loss regions extend radially from mid-radius to tip. 
Comparison of the relative flow measurements taken in the three- and four- 
bladed inducers provide the following conclusions: 

1) Measurements taken near the leading edge shows marked reduction 
in boundary layer growth, flow losses and radial inward velocity in the 
case of the three-bladed inducer. The losses near the tip are nearly 
halved from those of the four-bladed inducer and the "wake" type of 
profile observed near the mid-passage of the tip disappears in the case 
of the three-bladed inducer. 

2) Measurements taken near the trailing edge indicate that the 
relative velocity and pressure distributions (blade-to-blade) are 
similar for both of the inducers, even though there is appreciable 
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reduction in losses and relative flow diffusion in the case of the 

three -b laded configuration. 

* 

3) There is appreciable improvement in hub and wall static 
pressure distribution in the case of the three -bladed inducer. 

A method of measuring the three velocity components and the corres- 
ponding turbulent intensities has been developed for a stationary 
reference frame and has been reported in Refs. 12 and 20. This procedure, 
which utilizes three stationary hotwires located in the coordinate 
directions and located very near to the inducer trailing edge, has led 
to valuable information on the blade-to-blade variation of all the 
velocity and turbulence intensity components at that location. Major 
conclusions from these measurements are: 

1) The radial velocities obtained from the hotwire measurements 
are found to be of the same order of magnitude as the axial velocities 
throughout the flow passage. 

2) The blade -to -blade variation of axial, tangential and radial 
velocities measured at the exit by means of the hotwire probes are 
found to be nearly uniform. Similarly, the blade-to-blade variation of 
the relative velocities derived from the hotwire data is found to be 
nearly uniform. This is probably due to considerable wake diffusion 
that takes place between the trailing edge and the hotwire measuring 
station. 

3} The exit turbulence intensities are found to be rather uniform 
in the entire flow passage. The magnitudes of the turbulence intensities 
also reveal the highly turbulent nature of the flow in inducers. 

In an axial flow inducer, the interaction between the pressure 
surface and suction surface boundary layers result in an extremely 
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complex flow, especially near the outer half of the blade span. Since 
these interaction effects are very complex, experimental investigation 
assumes a very important role in this program. These effects are being 
studied at The Pennsylvania State University in a helical channel, 
specifically designed and fully instrumented for this purpose. The 
channel consists of constant -thickness blades of zero camber. It is 
proposed to carry out experiments at zero and other incidences. Details 
of the experimental program, the Reynolds equation and velocity profiles 
analysis valid for this flow and some preliminary measurements are given 
in Ref. 2. 


Methods and Means of Investigation 
As discussed in the previous sections, the three-dimensional viscid 
and inviscid effects of the inducer fluid flow makes meaningful 
predictions extremely difficult. An existing numerical procedure to 
simultaneously solve the three -dimensional equations of inviscid motion 
and continuity, developed by Cooper and Bosch (Ref. 4), is first utilized 
in an attempt to obtain the inviscid solution of the inducer flow field. 
The inducer geometry used for this numerical analysis is shown in Fig. 2. 
The incorporation of the dominant viscous terms into the equations of 
motion of Ref. 4 is then attempted. Use is made of empirically derived 
values of blade skin friction coefficient to derive the necessary 
viscous terms in the equation. The exact equations of motion, including 
all of the viscous terms, are extremely difficult to solve numerically. 

As an approximation, only the dominant viscous terms in the r-8-z 
directions are considered. 
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Developed view of 
sections 1, 4 


Figure 2. Inducer Geometry for Numerical 


Analysis . 
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Because of the iterative nature of the solution to the exact 
equations, a large amount of computer time and computer storage is 
usually required. In an attempt to reduce the solution convergence 
time, a method of initially determining input variables needed by the 
exact solution program is derived. 

Extensive measurements inside the blade passages are carried out 
not only to confirm the validity of the flow predictions, but also to 
gain a better understanding of the secondary motions of the flow. 
Extensive blade static pressure distributions are measured and compared 
with theory and a thorough survey of limiting streamline angles on the 
blade surface is also presented. A triaxial hotwire anemometer is used 
for measuring the three components of velocity, turbulence intensity and 
turbulence stress within the rotating inducer blade passage. The techni- 
ques used in this method are described in detail in a later section. 

The location of the experimental flow measuring stations are given in 
Fig. 3. Comparison of experimental results with results of the present 
theoretical analysis and results of previous investigations are dis- 
cussed. 

The three-b laded inducer was used in the experimental and theoreti- 
cal investigations of this thesis. 

Appreciation is expressed to E. Jordan and A. K. Anand for their 


assistance in instrumentation and measurements. 



Static pressure and limiting streamline 
angle measurements 

Flow survey station 

Figure 3. Location of Experimental Flow Measuring Stations 
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THEORETICAL ANALYSIS 


A thorough knowledge of all significant inviscid effects (blade 
blockage, flow turning, finite hub/tip ratio, etc.) and viscid effects 
(boundary layer growth, energy dissipation, etc.) is essential in the 
accurate prediction of the flow in all turbomachinery. Relevant to 
this, the availability of modem computers with large storage capacities 
and fast computation times greatly enhance the possibility of numerically 
solving the complete equations of motion. One of the early investi- 
gations in this area was made by Cooper and Bosch (Ref. 4) for the case 
of the three-dimensional inviscid flow through axial flow inducers. 
Application of this method of analysis to the Penn State inducer is 
given in this chapter. In addition, this chapter describes modifications 
to the Cooper-Bosch method which have been attempted to help reduce 
convergence time of the solution and provide a viscid solution capability 
based on empirically determined blade skin friction coefficients. A 
method of initializing the blade flow parameters as input to the Cooper- 
Bosch method has also been attempted in a search for a faster convergence 
to the solution. 


Exact Inviscid Analysis 

As mentioned above. Cooper and Bosch have developed a method of 
obtaining the exact inviscid solution of the inducer flow field. This 
three-dimensional analysis employs an iterative numerical procedure to 
solve the equations of motion expressed in finite -difference form. 
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General Equations and Method 

The nonlinear partial differential equations governing the flow in 
a rotating cylindrical coordinate system r, 0, z are: 


r momentum: 


g o 3 d „ 3U V 3U 3U 1 ( ~2 

“p37 +U 37 + F30 + W 37'F (V + rfl) 


+F = 0 
r 


( 2 ) 


0 momentum: !2.|| + u|^ + ^||-+w||- + ^ + 2Ufi + F Q = 0 


(3) 


z momentum: 


!o|E +U |W + V3W + w |W +f =0 
p 9z 9r r 90 9z z 


(4) 


continuity: 


U3U13V3W 
r 3r r 36 + 3z 


(5) 


Where W, V, U are relative velocities in the axial, tangential and 
radial directions respectively (Fig. 2) . F^, F^ and F^ are the 
components of the body forces including viscous terms, and are zero for 
the inviscid case considered in this section. In the Cooper-Bosch 
method, the above equations are rearranged to give residuals which are 
reduced to zero by a relaxation procedure. The total residual (RT) of 
one relaxation cycle is calculated by 

IMAX JMAX KMAX 9 9 9 9 

RT = ZEE [(RID + (R2 r + (R3 y + CR4)^]. JA (6) 

i=l j = l k=l 

where Rl, R2, R3 and R4 are the residuals calculated for the three 
momentum equations (2-4) and the continuity equation (5), and IMAX, 

JMAX and KMAX are the number of grid stations in the radial, tangential 
and axial directions which are used in the numerical analysis. 
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From this, the total RMS (root mean square) residual is defined as 


RMS = 


(RT) 2 

^ 4 • (IMAX) • (JMAX) • (KMAX) 


( 7 ) 


and is thus a measure of the degree of convergence between the iterated 
solution and the exact solution. 

The application of this method to the solution of the flow in the 
three-bladed Penn State inducer has been reported by Poncet and 
Lakshminarayana in Ref. 20. 

In applying this method to the Penn State inducer, the flow is 
assumed to be incompressible, and a grid of 7 x 7 x 26 is chosen to 
represent the blade passage. The flow geometry is shown in Fig. 2 , 

The boundary condition to be satisfied on the hub, annulus walls 

- V — f — y- 

and the blade surfaces is Q R ’ n = 0, where n is the direction normal 

-v 

to the channel boundaries and Q R is the total relative velocity. 

The first of the 26 axial stations corresponds to the upstream 
through-flow boundary where the initial conditions are applied. For 
the boundary value problem to be consistent, these initial upstream 
conditions must specify the three components of velocity and pressure, 
and the tangential velocity on the second axial station (which thus 
defines the swirl at the inlet of the inducer) . 

The last four axial stations correspond to the downstream flow- 
through boundary, and extend to about one-fifth of the chord length 

-v 

downstream of the trailing edge. With Q R - n = 0 to be satisfied on 
these stagnation stream surfaces, the set of boundary conditions for 
the problem is complete. 
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The Cooper-Bosch program has been suitably modified for use on 
the high speed digital computer system at The Pennsylvania State 
University. However, because of the iterative techniques employed in 
the Cooper-Bosch program, a large amount of computer time is usually 
required to converge to a satisfactory solution. For increased 
efficiency, the program has been compiled under a Fortran IV H level 
optimization procedure which reduces the time required for repetitive 
calculations, and production runs were submitted using the resulting 
object card deck. Output results from the computer program were placed 
onto 9 track, 1600 BPI (bits per inch) magnetic computer tape for future 
accessibility. 

Initial Input to the Exact Solution Program 

An investigation of available mathematical methods to solve the 
four simultaneous nonlinear partial differential equations governing the 
inducer flow revealed that there was no alternate method which would 
solve the equations more efficiently or effectively with a minimum of 
programming effort than the method described in Ref. 4. The next 
available approach for the speedier solution of the governing flow 
equations is the optimization of the input parameters of velocity and 
pressure which would allow faster convergence to the three-dimensional 
solution. Cooper and Bosch have derived an approximate solution in 
Ref. 4 to be used as an initial input to the exact program. This method 
derives the blade-to-blade average quantities using axisymmetric 
equations, then uses these quantities in a blade-to-blade solution of an 
integrated form of the scalar momentum equation in the tangential 
direction. The flow parameters derived by this method were used in 
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Ref. 20. However, an alternative method of developing the initial input 
flow parameters has been attempted in the present analysis, and is 
discussed in detail in the following two sections. 

Douglas -Neumann Analysis . - The initial estimation for the velocity 
and static pressure distribution throughout the inducer flow passage is 
calculated by the two-dimensional Douglas -Neumann program described in 
detail in Ref. 5. The technique employed by the Neumann program to 
solve a particular fluid flow problem is to use source distributions of 
appropriate strength on the surface of the blade profile in such a way 
that the flow normal to the surface of the body is either zero or 
prescribed. When the Neumann boundary condition is applied, an 
integral equation in source strength £ is obtained 

-C • n = £(s) + / C(x)A(x,s) dx (8) 

body 

— >• — ■ 4 ' 

where A(x,s) = n • C(x,s) and C is the onset flow. C(x,s) is the 
velocity at a surface point s due to a unit source at x. The solution 
for the general case of a lifting cascade at any angle of attack is 
calculated by superposition of three ,T basic flows” in such a way that 
the correct angle of attack is obtained and the Kutta condition is 
satisfied. The ”basic” flows are: flow at zero angle of attack, flow 

at 90° angle of attack, and circulatory flow for each cascade. Super- 
position of solutions is possible because the potential equation is 
linear and the boundary condition on the cascade blade is homogeneous. 
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In the Douglas -Neumann results, velocities and static pressure 
coefficients are normalized with the modulus of the average onset flow 
velocity 

C. n + C 

inlet exit rCk ^ 

2 ’ ) 

The listing of the Douglas -Neumann program is given in Ref. 5. 

Quasi-Three-Dimensional Modification. - The investigation of the 
flow around an isolated airfoil in a contracting or diverging stream is 
presented in Ref. 13. This analysis provides a simple method of modi- 
fying the two-dimensional Douglas -Neumann flow solutions to account for 
the three-dimensional effect of the converging or diverging streamlines. 
An expression for static pressure coefficient on the airfoil surface is 
derived as a function of channel slope, two-dimensional static pressure 
coefficient, and the Fourier coefficients of the blade profile. The 
analysis utilizes thin airfoil theory approximations and assumes that 
thickness effects are the same as in plane flow. The mean flow is 
assumed to be inviscid, steady, and incompressible, and the variation 
of channel height is assumed to vary linearly from leading to trailing 
edge while the length of the contracting section is assumed to be the 
same as the axial projection of the blade. In the present application, 
the expression for static pressure coefficient has been modified in an 
attempt to represent the flow about a row of two-dimensional infinite 
cascades. TTiis quasi -three-dimensional approach has been applied to 
the two-dimensional results obtained from the Douglas -Neumann analysis 
for the Penn State inducer. The effect of the converging channel as 
determined by the above analysis on the Neumann solution for the Penn 
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State inducer is essentially to decrease the blade static pressure 
near the trailing edge. 

A comparison of the radial variation of axial and tangential 
velocities calculated by the above method with the experimental results 
of Ref. 20 shows close agreement (Fig. 4) . The agreement between the 
measured tangential velocity and inviscid prediction may be fortuitious, 
since the axial velocity predicted at the same location is considerably 
different from the measured values. 

Using the input parameters of velocity and pressure derived from 
the preceeding analysis results in a lower total RMS (root mean square) 
residual than with the previous method of initializing the input 
variables. As an example, the final RMS residual for the inviscid 
results of Ref. 20 was 0.12450 after 68 relaxation cycles, whereas a 
similar value is obtained using the present analysis in 10 relaxation 
cycles. This amounts to a considerable saving in computer time. 

Twenty- five iteration cycles has reduced the RMS residual to 0.10579, 
indicating that a faster convergence to the solution should be possible. 
Further investigation should be carried out to confirm the effectiveness 
of the input analysis as an alternative to the Cooper-Bosch approximate 
solution method. 

In a further attempt to decrease the convergence time, the exit 
flow angle was allowed to change depending upon the tangential and axial 
velocities calculated at the inducer trailing edge. Since the exact 
downstream boundaries are not known in this type of problem, it was 
hoped that by allowing the downstream boundaries to adjust themselves 
and thereby unload the blade trailing edge, a more exact definition of 
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the downstream streamlines would result in lower RMS residuals. Cooper 
and Bosch suggest a similar technique as a means of reducing RMS 
residuals in their recommendations for future work. 

Since the extension of the stagnation stream surfaces downstream 
have, been constructed to be uniformly periodic with a spacing of 2 tt/N 
(N being the number of blades) , the values of velocity and pressure at 
the downstream tangential channel boundaries should be equal. This 
condition is applied at the blade trailing edge after each iteration 
cycle. If the pressure and suction surface parameters differ with each 
other at the trailing edge grid point, the average value is used in the 
residual calculations. If the axial and/or tangential velocities at 
the trailing edge diverge significantly from the design values during 
the iteration process, then the flow exit angle, defined by $ = tan 
at the trailing edge, is recalculated and is used to redefine the down- 
stream stagnation stream surfaces. This method also has the advantage 
of automatically forcing the Kutta-Joukowski condition for the blade 
pressure distribution to be satisfied. Changes made to the original 
Cooper-Bosch program can be seen in Appendix B and are concentrated in 
subroutine "MAIN". A flow chart diagram of the Cooper-Bosch program, 
including the modification discussed above, is given in Appendix A. 

Viscid Analysis 

In addition to the attempts to improve the convergence of the 
exact inviscid solution, a method of incorporating viscid effects into 
the governing equations of motion has also been investigated. 


=l< 



General Equations and Theoretical Method 
The viscid equations of motion are: 
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Molecular viscosity terms have been neglected in these equations. 
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Comparing these equations with the momentum equations used by Cooper 
and Bosch, the following expressions for FR, FT, and FZ (the exact 
program variables for viscous loss terms) can be given as: 
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Since the stagger angle is very large, these viscous terms can be 
approximated by retaining the dominant terms as well as neglecting the 
normal shear stresses, resulting in: 
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(17) 

(18) 
(19) 


The distribution of shear stress is assumed to be linear across 
the flow passage from pressure surface to suction surface. The values 
of wall shear stresses are assumed to be known from previous experi- 
mentation. Skin friction coefficient C^ for a four-blade flat plate 
helical channel is given in Ref. 2. The results, summarized in Fig. 5, 
are considered to be valid for the three-blade inducer under considera- 
tion. Interpolation of the curves in Fig. 5 for a given blade surface 
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Figure 5. Skin Friction Coefficient vs . Reynolds Number 
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grid location under consideration gives a value of wall shear stress 
T w = • j- p Q R for the appropriate Reynolds number R e = Q r t/\> of the 

flow at that point, where Q” r is the average relative velocity across 
the flow passage as derived by the Cooper-Bosch relaxation procedure. 
Calculation of x at each grid location within the flow passage allows 
the derivatives of equations 17-19 to be calculated by finite-difference 
methods . 

An additional requirement placed on the viscid analysis is to 
satisfy the viscid boundary condition which requires that all components 
of velocity are zero at the blade surface. 

The changes to the original Cooper and Bosch exact program 
necessitated by the inclusion of the viscous loss terms are made in 
subroutines "MAIN", "DLOSS" and "RESID". Flow chart diagrams for the 
modified subroutines of the Cooper-Bosch program are given in Appendix 
A. A complete Fortran listing of the modified Cooper-Bosch program is 
given in Appendix B. 

Input and Solution 

The input variables and formats for the modified viscid analysis 
program are identical to the original Cooper-Bosch program, with the 
exception of including a set of curves to define blade skin friction 
coefficient (C^) vs. Reynolds number (R g ) for various reference 
tangential locations throughout the inducer channel. Input values are 
taken from log-log plots similar to that shown in Fig. 5. Straight line 
approximations for the reference data are required. A definition for 
each of the additional input quantities follows. For the exact format 
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in which these parameters must be coded, consult the program Fortran 

listing in Appendix B. 

NCURVE: Number of vs. R g reference curves used as input. Each 

curve must correspond to a specific blade tangential 
location. Maximum of 3, minimum of 2. 

THETA(I) : Tangential location, in degrees from leading edge, where 

a specific C_ vs. R curve applies. The array index (I) 
increases from 1 to NCURVE, proceeding from leading edge 
to trailing edge. 

REREF1 (I) , 

REREF2(I): Minimum and maximum values respectively of R used in 

straight line approximation of vs. R g curve. 

CFREFP(1,I) , 

CFREFP (2 , I) : values corresponding to REREF1 (I) and REREF2 (I) 

respectively on the blade pressure surface. 

CFREFSfl, I) , 

CFREFS(2,I): values corresponding to REREFl(I) and REREF2(I) 

respectively on the blade suction surface. 
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The viscid modifications which have been discussed are activated 
in the computer program when the appropriate value of fluid kinematic 
viscosity is used as an input parameter. If zero viscosity is coded, the 
modified program will automatically revert to an inviscid analysis as 
represented in the original Cooper-Bosch program. Preliminary running 
of the modified viscid program indicates an increase in computer time of 
approximately two to three times more than a corresponding inviscid 
analysis run. This increase in computation time is due to the calculation 
of the viscous loss terms FR, FT and FZ at each grid point location 
throughout the duration of one relaxation cycle, which may involve 
several thousand iterations of the flow parameters in order to reduce 
the RMS residual from the previous cycle. The efficiency of the computer 
programming can certainly be improved upon in future study. 

Comparison of the results of the Cooper-Bosch inviscid and viscid 
analyses and their relationship to experimental measurements are 
discussed in a later chapter. 
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EXPERIMENTAL EQUIPMENT, METHODS AND TECHNIQUES 

The primary goal of performing the following experimental program 
was to investigate the flow characteristics such as velocity, turbulence 
characteristics and static pressure of the relative flow inside a three- 
bladed inducer. The importance of this experimental data for a better 
understanding and prediction of the flow has been mentioned in the 
introduction . 

A secondary goal was the determination of the effectiveness of the 
triaxial rotating hotwire as a method of measuring mean and fluctuating 
velocities and turbulence stresses within the inducer blade passage. 

Three -Bladed Inducer 

The experimental investigation was performed on a three -foot 
diameter axial flow inducer with three equally spaced blades. The test 
facility is pictured in Fig. la. Design of the blades is by the mean 
streamline method of Wislicenus (Ref. 25) . The inducer was operated at 
450 rpm, which was determined to an accuracy of 0.1 rpm by means of a 
photocell circuit with rotating calibrated disk and displayed on an 
electronic counter. Important parameters of the inducer are' as follows: 


Number of Blades 

3 

Hub/Tip Ratio at Outlet 

0.50 

Hub/Tip Ratio at Inlet 

0.25 

Radial Clearance 

0.0625" 

Inlet Flow Coefficient (Design) 

0.065 

Blade Chord at r/r t = 1.0 

82.96" 

Blade Chord at r/r t = 0.75 

63.18" 
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Blade Chord at r/r t = 0.50 

49.94" 

Solidity at r/r t = 1.0 

2.15 

Solidity at r/r t = 0.75 

2.21 

Solidity at r/r^ = 0.50 

2.61 

Reynolds Number Based on Tip Radius 

7.0 x 10 


The design values of blade and flow angles at inducer inlet and exit 
are given in Fig. lb. 

The use of the three-bladed inducer for the continued experimental 
investigation defined in this report is a result of conclusions reached 
by prior investigations described in Ref. 20; namely, it has appreciably 
better performance than a similar four-bladed inducer tested at the same 
flow coefficient. 

Blade static pressure measurements were obtained with the use of 
hypodermic steel tubing of .063" ID imbedded in the blade at ten 
separate pressure and suction surface locations. The pressure measure- 
ments at each location were carried out at five radial stations by 
utilizing .063" diameter taps drilled at equally -spaced intervals from 
tip to hub. The approximate radial locations of the pressure taps are 
shown in Fig. 3 and the actual experimental radial and chordwise locations 
of the pressure taps are given in Table 1 and Table 2 respectively. 

Velocity and turbulence measurement stations within the blade 
passage have been previously used (Ref. 20) and were constructed by 
cutting tangential slots in the hub wall at the locations shown in 
Fig. 3. 

The inducer was statically and dynamically balanced at facilities 
in the Garfield Thomas Water Tunnel of The Pennsylvania State 
University. 



Table 1. Radial Location of Blade Static 


Pressure Measurement Stations* 


Tangential 

Measurement 

Station 

1 

Radial Measurement Station 
2 3 4 

5 

1 

.291 

.466 

.644 

.815 

.985 

2 

.295 

.466 

.637 

.808 

.985 

3 

.322 

.479 

.654 

.819 

.985 

CD 4 

L) ^ 

.342 

.500 

.664 

.823 

.985 

aS 

m 

u S 

.356 

! .514 

.671 

.823 

.985 

in 

2 6 

.370 

.521 

.678 

.829 

.985 

3 

t/> 

</> 7 

<u ' 

.390 

.541 

.689 

.835 

.985 

U 

Du 

8 

.411 

.555 

.699 

.842 

.985 

9 

.438 

.575 

.719 

.849 

.985 

10 

.473 

.609 

.726 

.863 

.985 

1 

.288 

.466 

.640 

.815 

.985 

2 

.301 

.479 | 

.644 

.823 

.985 

3 

.322 

.486 | 

.658 

.823 

.985 

4 

.342 

.500 , 

.664 

.823 j 

.985 

<u 

o 

4-1 D 

.356 

.514 

.671 

.829 

.985 

U 

3 

w 6 

.370 

.527 

.678 

.835 

.985 

§ 

a 7 

.390 

.541 

.685 

.835 

.985 

g 

w 8 

.397 

.548 

.692 

.835 

.985 

9 

.425 

.555 

.699 

.835 

.980 

10 

.435 

.561 

.692 

.842 

.983 


*Radial locations nondimensionalized as r/r^. 
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Table 2. Chordwise Location of Blade Static 
Pressure Measurement Stations* 


Tangential 

Measurement 

Station 

1 

Radial Measurement Station 
2 3 4 

5 

1 

11.3 

8.9 

6.5 

5.7 

4.9 

2 

15.5 

14.3 

13.0 

12.0 

11.0 

3 

24.0 

24.0 

24.0 

23.8 

23.5 

<D 4 

o 

34.0 

36.0 

38.0 

38.0 

38.0 

aJ 






D 5 

44.0 

48.0 

52.0 

52.3 

52.5 

CO 






g 6 

51.0 

55.5 

60.0 


61.0 

D 






2» 7 

CD 

58.0 

63.0 

68.0 

69.0 

70.0 







Cl, 






8 

66.0 

72.0 

78.0 

78.5 

79.0 

9 

74.0 

i 

80.0 

86.0 

86.5 

87.0 

10 

82.0 

88.0 i 

94.0 

95.0 

96.0 

1 

15.0 

13.3 

11.5 

11.0 

10.4 

2 

19.5 

18.9 

18.2 

| 

17.6 

17.0 

3 

28.4 

29.2 

30.0 

i 30.0 

30.0 

4 

45.0 

47.5 

50.0 

50.5 

51.0 

a) 






CJ 






rf 5 

q-l ° 

54.0 

57.0 

60.0 

61.0 

62.0 

fn 






d 






co 6 

62.0 

66.0 

1 70.0 

70.8 

71.5 

jd 






o 






■H 7 

+J 7 

69.0 

73.5 

78.0 

78.5 

79.0 

o 






co 8 

75.0 

80.5 

86.0 

86.0 

86.0 

9 

81.0 

86.5 

92.0 

92.3 

92.5 

10 

85.0 

91.0 

97.0 

97.0 

97.0 


*Chordwise locations are expressed as percent chord from blade leading 
edge . 
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Blade Static Pressure Measurement 
An extensive experimental investigation of the blade static 
pressure distribution has been undertaken to help provide a check on 
prior theoretical analyses and useful information for future theoretical 
development. 

Equipment Used 

A schematic diagram of the blade static pressure test setup is 
shown in Fig. 6. The equipment used to measure the blade static 
pressure distribution of the three-bladed inducer is as follows: 

Scanivalve . - The scanivalve, a scanning type pressure sampling 
valve for measuring multiple pressures, was mounted in the rotating hub 
section of the inducer. The scanivalve incorporates a fluid wafer 
switch for time-sharing one pressure lead with up to twenty-four (24) 
unknown pressures, and is stepped by a rachet-geared solenoid. A 
solenoid controller used push button pulse length feedback and increased 
drawing voltage to step the solenoid driven scanivalve. The controller 
was equipped with a 24-division indicator dial which allowed monitoring 
of the static pressure station under consideration. 

Three Channel Pressure Transfer Device (PTD) . - A 3-channel pressure 
transfer device was used to transfer the static pressure measurements 
from the rotating reference frame of the three-bladed inducer to the 
stationary reference frame. Each channel was made airtight by the use 
of double -sealed ball bearings, and pressure leakage was prevented by 
use of O-rings and plastic sealers. The PTD was mounted on a stand 
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Figure 6. Schematic Diagram of Blade Static Pressure Test Setup 
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outside the rotor assembly and was housed in a streamlined cowling to 
reduce any interference on the incoming flow. 

Slip-Ring Unit . - An eight -channel slip-ring unit was used to 
conduct electrical signals from the stationary reference frame to the 
rotating reference frame of the scanivalve. Electrical continuity was 
provided by carbon brushes in contact with a rotating commutator aligned 
along the inducer's rotational axis. The slip-ring unit was mounted on 
the pressure transfer device, and all electrical and pressure connections 
were transferred through a hollow shaft and flexible couplings to the 
nose cone of the inducer. 

Peripheral Equipment. - A transistorized 30 volt D.C. power supply 
was used to provide voltage to the scanivalve and solenoid controller 
units. A micromanometer graduated in 0.001" divisions was used to 
measure the blade static pressure. 

Procedures and Techniques 

Flexible vinyl tubing of 0.063" inside diameter was used to connect 
the ten suction and ten pressure surface stations to the available tubes 
on the scanivalve. The vinyl tubing was also used to connect the 
collection tubulation of the scanivalve to the measurement channel of 
the pressure transfer device. The manometer was similarly joined to this 
channel. Electrical connections from the scanivalve were transferred 
through the slip-ring unit to the solenoid controller and power supply. 
The blade static pressure measurement test setup is shown in Fig. 7. 
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Masking tape was used to cover the radial pressure taps not under 
consideration in a particular run. 

Once the inducer speed was adjusted to 450 rpm, the static pressures 
associated with the ten pressure and ten suction surface stations at a 
particular radial passage location could be measured by stepping the 
scanivalve through its range of operation. The dial on the solenoid 
controller would indicate which station pressure was being measured. 

Each step provided the blade static pressure of a different chordwise 
station. Pressure readings were measured on the micromanometer to an 
accuracy of 0.001". 

Since the blade static pressure measurements were taken on the 

rotating blade, it was necessary to apply a centrifugal force correction 

to obtain the static head. If h is the height of the water column 

m 

measured by the manometer, the actual static head is given by 


P r, 2 

m , ft , ^ 
— h + ^ — (r 
P m 2g Q 


r o 2) 


( 20 ) 


where p m is the density of the manometer liquid, r Q is the radius of the 
rotating shaft used in the pressure transfer device, and r is the radius 
of the static pressure tap under consideration. 

From this, the blade static pressure coefficient is defined by 


= 

r s 


2g o h 

U. 2 


( 21 ) 


and is calculated for all pressure measurement stations. 
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Blade Limiting Streamline Angle Measurement 
The measurement of the blade limiting streamline angle (a) , which 
is the limiting position of the streamline as the blade surface is 
approached, is a further attempt to define the flow phenomena within the 
blade boundary layer, including nature and magnitude of the radial flows 
and the direction of the wall shear stress. The information gained will 
help establish the extent of three-dimensionality in the inducer flow 
and will be valuable in developing a viscid theory for the prediction of 
rotating boundary layer characteristics. 

Equipment, Procedures and Techniques 

The blade measurement stations are identical to those used for the 
blade static pressure measurements of the previous section. The method 
and equipment for measurement are essentially the same as that used in 
Ref. 10. An ammonia transfer device (ATD) was placed inside the hub 
section along the axis of rotation. For a specific blade measurement 
station, the .063” flexible vinyl tubing associated with that location 
was attached to the ATD. Likewise, a thin strip of ozalid paper was 
attached radially on the blade surface adjacent to the measurement 
station. 

Once the inducer had been rotated to the required 450 rpm, a small 
amount of ammonia gas was injected at low mass flow rate into the ATD 
and allowed to penetrate through the radial taps on the blade. The 
resulting ammonia trace on the ozalid paper was then measured to deter- 
mine the limiting streamline angle. This procedure was repeated for all 


measurement stations. 
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Three -Sensor Rotating Hotwire Measurement 
Initial feasibility investigation into the use of the hotwire 
anemometer in the rotating flow passage of the three -blade Penn State 
inducer has been mentioned in Ref. 2. One problem experienced had been 
the inability of the relatively crude slip-ring unit to maintain the 
continuity of the hotwire circuitry. It did not allow suitable operation 
of the hotwire for extended periods of time to allow meaningful measure- 
ments to be obtained. The present study was an attempt to improve upon 
the accuracy and longevity of the initial investigation and prove the 
suitability of hotwires in the measurement of the relative mean and 
fluctuating velocities in a rotating environment. 

Equipment Used 

A schematic diagram of the rotating hotwire test setup is given 
in Fig. 8. A detailed description of the equipment used in the measure- 
ment of the relative mean and fluctuating velocities within the rotating 
passage of the three-bladed inducer follows: 

Triple -Sensor Hotwire Probe. - A subminiature triaxial probe 
designed for boundary layer flows was used in the experimentation (Fig. 

9) . The wire is 3ym diameter copper plated tungsten with a length/ 
diameter ratio of approximately 300. The probe was attached to a 
specifically designed probe support for use in traversing the inducer 
flow passage (Fig. 10a). 
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Figure 8. Schematic Diagram of Rotating Hotwire Test Setup 


















Figure 10a. Triaxial Hotwire Probe 
and Support 


Figure 10b. Triaxial Probe 
Mounted in Inducer Hub 
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Three-Channel Constant Temperature Hotwire Anemometer .- Two dual- 
channel constant temperature hotwire anemometers were used to provide 
the three-channel capability necessary for these measurements. The 
original anemometer circuitry has been given in Ref. 20 and is not 
reproduced here. 

Mercury Slip-Ring Unit .- A ten-channel mercury slip-ring unit was 
utilized in transmitting the hotwire signals from the rotating reference 
frame of the inducer to the stationary hotwire anemometers. The slip- 
ring unit exhibits the smallest and most stable resistance in the 
transfer of measuring signals from the rotating electrical elements to 
the stationary electrical conductors. Contact between the rotating 
wires and the stationary contact screws is made through a round contact 
disc to which the rotating wire is connected, rotating in mercury. 
Triple-distilled mercury was used to provide the greatest conductivity 
and the lowest noise level distortion possible. 

Peripheral Equipment .- The A.C. voltage signals obtained from the 
hotwire anemometer were visualized on a four-channel storage oscilloscope. 
Instantaneous mean velocity D.C. voltage readings from the anemometers 
were displayed on a digital voltmeter. 

The fluctuating voltage (A.C.) signals were processed through a 5.0 
KHZ low-pass filter driven by a 15-volt regulated power supply, which 
was used to cut off the high frequency noise which may have entered the 
circuit . 

Mean-square values of the A.C. voltages were obtained by passing 
the signals through a true RMS voltmeter and subsequently through a 
manually controlled signal integrator. The mean-squared voltage was 
displayed on a digital voltmeter. 
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A sum-and-difference circuit was utilized to obtain the sum and 
difference between the three hotwire signals needed for the turbulence 
intensity calculations. 

A signal generator was used for sinewave generation to determine 
gains throughout the hotwire circuitry and the accuracy of the 
associated peripheral equipment. 

Calibration Equipment . - A low-turbulence calibration tunnel was 
used for the hotwire calibration. The horizontal wind tunnel has a 
test cross-section of 1-1/2" x 1-1/2" and operates within the range of 
air velocities of 0 to 300 feet per second. The calibration velocities 
were measured with a nonshielded pitot tube and the micromanometer 
described previously. 

Procedures and Techniques 

Measurements were taken at two axial stations, corresponding to 
approximately 33% and 90% of the blade chord (Fig. 3) . Various 
velocity measurements have previously been performed at these stations 
(Ref. 20) and thus a comparison of hotwire experimental results with 
these prior investigations are possible. 

Six radial stations (corresponding to r/r^ values of .973, .945, 
.890, .781, .671, .548) at station 1 and five radial stations (r/r^ 
locations of .973, .945, .890, .781, .671) at station 2 were traversed 
at several tangential intervals within the blade passage in an attempt 
to get an accurate and detailed appraisal of the flow velocities. 
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turbulence intensities and stresses in these regions. It was not 
possible, however, to obtain measurements extremely close to the blade 
surfaces due to the limitations caused by the blade curvature. 

The three-sensor probe was attached to a ten-inch section of 
aluminum tubing and fixed in a particular radial and tangential location 
by a coupling mounted in the inducer hub (Fig. 10b) . The probe was 
accurately aligned in the tangential direction with the aid of the guide 
vane attached to the probe's adjustable protection pin. Orientation of 
the three individual hotwires was measured with respect to the (R, 0, Z) 
coordinate system (Fig, 9) by utilizing a linearly -calibrated scale eye- 
piece in a 30-X microscope. The direction cosines of this orientation 
were then calculated, as were direction cosines of the two arbitrary 
normals to each wire (Table 3) . These values were used in the governing 
hotwire equations derived in the next chapter. 

The experimental setup for the hotwire measurements is shown in 
Fig. 11. With the probe in position, the inducer was started and 
rotated to 450 rpm. The corresponding mean D.C. voltages of the three 
hotwire channels E E^ were recorded, in addition to the statistical 
properties of the fluctuating voltages 

~2 ~ 2 ~2 T 77 7 TT 

e l « e 2 ’ e 3 » + e 2 ) * < e i ' e 2 ) » 

(6 1 + V 2 ’ (S r e 3 )2 ’ (e 2 + e 3 )2 and (e 2 " e 3 )2 ‘ 

The time -averaged voltages were obtained over an integration of 100 
seconds. The inducer was then stopped, the probe was moved to another 
location, and the procedure repeated until the flow field was entirely 




Table 3. Direction Cosines Used in Hotwire Analysis 


00 


Angle 

Orientation 

Hotwire 1 

Hotwire 2 

Hotwire 3 

(Fig. 9 and 
Eqns. 22-24) 

Angle 

(degrees) 

Direction 

Cosine 

Angle 

(degrees) 

Direction 

Cosine 

Angle 

(degrees) 

Direction 

Cosine 

ax-6 

a i 

117° 24.5' 

-.46034 

120° 59.2' 

-.51486 

130° 

26.5' 

-.64865 

AX-R 

b i 

128° 53' 

-.62773 

49° 56.5' 

.64357 

94° 

39' 

-.08108 

AX-Z 

c i 

51° 07' 

.62773 

55° 30.5' 

.56634 

139° 

21.5* 

-.75676 

Nl-0 

a 2 

90° 

0 

90° 

0 

90° 

0 

Nl-R 

b 2 

45° 

.70711 

48° 39.1' 

.66063 

173° 

53' 

-.99431 

Nl-Z 

C 2 

45° 

.70711 

138° 39.1' 

-.75071 

83° 

53' 

.10653 

N2-6 

a 3 

152° 35.5' 

-.88775 

149° 00.8' 

-.85728 

139° 

33.5' 

-.76109 

N2-R 

b 3 

71° 00.2' 

.32551 

112° 44.5' 

-.38651 

86° 

02.3' 

.06910 

N2-Z 

C 3 

108° 59.8' 

-.32551 

109° 53.1' 

-.34013 

49° 

50.5' 

.64496 
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surveyed. Station 1 and Station 2 were surveyed similarly, expect that 
no turbulence measurements were recorded at station 1. 

The resultant voltage measurements from the three- channel rotating 
hotwire experiment were converted to mean velocities U, V, W and 

~ T ~ 

turbulence quantities u , v , w , uv, uw, vw from the appropriate cali- 
bration curves and the applicable equations derived from the analysis 
of the next chapter. 

The data reduction was accomplished in a computer program written 
to solve, for all flow stations considered, the resulting three 
simultaneous mean velocity equations and six simultaneous turbulence 
velocity equations. The high speed digital computer at the Penn State 
Computation Center was used in this task. 
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DERIVATION OF THE HOTWIRE EQUATIONS 

The derivation of the equations for the three sens or -hotwire is a 
completely general method and is based on the hotwire configuration 
shown in Fig. 9. 


Direction Cosine Method for Nonorthogonal Probes 
Consider one hotwire sensor with respect to the (R, 0, Z) coordinate 
system. Hotwire (1) has an orthogonal coordinate system (AX, Nl, N2) 
associated with its orientation. This (AX, Nl, N2) coordinate system 
can be transformed to the (R, 0, Z) coordinate system by: 

(Q + q) ax = a x (V + v) + ^ (U + u) + c 1 (W + w) (22) 

(Q + q) nl = a 2 (V + v) + b 2 (U + u) + c 2 (W + w) (23) 

(Q + q) n2 = a 3 (V + v) + b 3 (U + u) + c 3 (W + w) (24) 

where Q ax > Q n ^ and Q n2 are the mean velocities associated with the 
(AX, Nl, N2) coordinate system and U, V and W are the transformed mean 
velocities in the (R, 0, Z) coordinate system, q, u, v, w are the 
fluctuating components. The coefficients a b^, c^, etc. are the 
applicable direction cosines between (AX, Nl, N2) and (R, 0, Z) . The 
specific direction cosines for the triaxial hotwire probe used in the 
experimental studies of this thesis are given in Table 3. 

The effective cooling velocity sensed by hotwire (1) is known 


(Ref. 22) to be: 
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CQ + q) x - ICQ + <oh * (Q + q)^ 2 + K x 2 (Q + q)^] 172 (25) 

where is the associated correction factor for deviation from the 
cosine law. 

Substituting equations 22-24 into equation 25 and expanding the 
resulting expression to reflect the instantaneous velocity components 
as the sum of the mean (Q^, U, V, W) and fluctuating (q^, u, v, w) 
velocity components, we get: 

Qj_ + q x = ta 4 (V 2 + v 2 + 2Vv) + b 4 (U 2 + u 2 + 2Uu) 

+ c 4 (W 2 +w 2 + 2Ww) + d 4 (UV + uv + Uv + Vu) 

+ e^(VW + vw + Vw + Wv) 

+ f 4 (UW + uw + Uw + Wu)] 1/2 (26) 


where the constant coefficients are defined as 

2 2 „ 2 2 
a 4 " a 2 * a 3 * K 1 a l 


. .2 .2 „ 2 . 2 

b 4 “ b 2 + b 3 * K 1 b l 


2 2 2 2 
C 4 * C 2 * C 3 * K 1 C 1 


d 4 = 2 (a 2 b 2 * a 5 b 3 * K 1 .jbp 
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e 4 = 2 (a 2 C 2 * a 3 C 3 * K 1 Vl> 


f 4 - 2 (b 2 c 2 . b 3 c 3 4 K 1 ‘b 1 c 1 ) 


The right side of equation 26 can be linearized by 
(1 + A) 1/2 = 1 + A/2 - A 2 /8 + . . . 


(27) 


where A is considered small. Thus, equation 26 can be approximated to 
the following expression: 

Ql Mj s a 5 (V + v) + a fi (U + u) + a ? (W + w) 


,,2 2 „,2 2 11M 

r U ^ u . , .W w . ^ ,UW uw . 

*8^ V ^ a 9^v + V ^ a !0 C V + V ^ 


Uu Ww .Uw Wu. ^ n . . 

+ a ll- + a 12 V” + a 13 + V" 5 + 0(e) 


(28) 


where the constant coefficients are defined by 


a s = ’ r ~ a 7 


a 6 = d 4 /2 


a 7 = e 4 /2 / a 4 


a 8 = ' rT 4 




a 9 = /_a 7 


[!i_. i ( !i) 2 ) 

4a,, 8 '■a. J J 

4 4 


a io = ^ 


f. , d e 

[ 2a7 " T ( a7^ ( i7 )] 
4 4 4 
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a ll = /I T 


[— - — (— ) 2 ] 
*4 4 J 


a 12 = 


£ - 4 ^ ^ 


a 13 = a 10 


Taking the time -average of equation 28, we get: 


U 2 2 

Q 1 = a 5 V + a 6 U + a 7 W + a 8 C V~ + T* 


,W 2 w 2 . f UW . uw,. 

+ a 9*-V + V -* + a 10 *-V V ^ 


(29) 


2 2 — 

If we consider the fluctuating quantities u , w and uw to be small, 
then the equation for the mean velocity sensed by the hotwire can be 
approximated by: 


a 5 V * a 6 U * a 7 W * a 8 T 


w uw 

+ a 9 V + a l0 V 


(30) 


The difference between equations 28 and 29 is the fluctuating 
velocity sensed by the hotwire: 


u 2 u 2 

*1 = a 5 V + a 6 U + a 7 W + a 8 3 


2 2 — 

.w w . .uw uw. . Uu 

+ a 9 *-V " V ^ + a i0 ( V " V ^ a ll V 


Ww , ,Uw Wu, 

+ a l2 V a i3 + V -* 


( 31 ) 



Expressions similar to equations 30 and 31 for Q 2 » Q 3> q 2 and q 3 
corresponding to the mean and fluctuating velocities sensed by hotwires 
(2) and (3) can be derived. 


Relating Voltages and Velocities 
Application of King’s Law (Ref. 6) for hotwire (1) relates 
instantaneous velocity (= + q^ and instantaneous voltage (= E^ + e 1 

by the expression 

(Ej + e^ 2 = E q 2 + E 1 ✓ Q x + q x (32) 

where is the slope of the hotwire calibration curve, E q is the hot- 
wire voltage at zero velocity, E^ is the mean voltage and e^ is the 
fluctuating component. 

Expanding and linearizing equation 32 for small q^/Q^, we get 
E x 2 + 2E l6l + e x 2 = E q 2 + B 1 /Qj' (1 + \ q 1 /Q 1 - ... ) (33) 


The time -average of equation 33 relates the mean velocities and 

2 

voltages. Neglecting e^ , we find 


E 1 2 • E o 2 * B 1 ^ 


( 34 ) 


Subtracting equation 34 from equation 33 and neglecting e.. , 


then 
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or 



( 36 ) 


Similar expressions can be derived for hotwires (2) and (3) . Rewriting 
equation 36 as q^ = L^e^, the time-averaged fluctuating quantities can 
be found to be 


< 2 t 2 2 > 
“ii ' L i e i 


q i q 2 L 1 L 2 ®l e 2 

, 2 2 2 



q 2 " L 2 ®2 

> and < 

V 3 * L 1 L 3 e l e 3 [ 

2 T 2 2 



q 3 " L 3 ®3 J 


, q 2 q 3 = L 2 L 3 e 2 e 3 / 


2 2 2 

The values of e^ , e^ and e^ are obtained directly from the experi- 


mental hotwire measurements. The values of e.e-, e.,e„ and e 0 e^ are 

1 2 * 1 3 2 3 


derived as follows: 


V 2 T [te l * e 2 )2 

■ (6 1 

- e 2> 2 ! 

e l e 3 T lte i * V 2 

“ (6 1 

- 

e 2 e 3 T Ue 2 ' e 3 )2 

- ( g 2 

- e 3 ) 2 ] 


(38) 


The mean -squared voltages within the brackets are obtained by 
utilizing the sum-and -difference circuit in the experimental hotwire 
measurement (Fig. 8) . 
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Mean Velocity Calculation 

Equation 30 for hotwire sensor (1) and similar equations derived 
for sensors (2) and (3) form a set of three nonlinear simultaneous 
equations in three unknowns: 


2 2 

U W UW 

Ql = a 5 V + a 6 U + a ? W + a g ^ + a g ^ + a 1Q — 


Q 2 = b 5 V ♦ b 6 U * b 7 W ♦ b 8 £ * b 9 ”1 ♦ b 10 f- 


- „ „ w U 2 W 2 UW 

% = c V * c 6 U * C 7 W * c 8 r * c ♦ c 10 5- 


where the coefficients a^, b 5 , etc. reflect the appropriate combination 
of direction cosines, K factor, etc. 

Values of C^, Q 2 and are known from application of equation 34 
to the D.C . hotwire voltages obtained from the experimental hotwire 
measurements. Solution of equations 39-41 thus give U, V and W for each 
measurement location considered. 

The Newton -Raphson method provides an iterative procedure for 
solving a nonlinear system of equations involving n real functions and 
n real variables. Details of the method are given in Ref. 3. Applying 
this method to the three simultaneous nonlinear equations 39-41, we can 
define the expression 



(42) 

(43) 
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f 3 ( u > V * W) 


c 5 v 


C 6 U 


C 7 W 


+ c r 


l LT 
V 


+ c, 


9 V 


+ c 


10 


™ - Q 
V ^3 


(44) 


After determining the partial derivative expressions 


3^ 3^ 

"3U ’ “3V ' “§W 

' i = 1,3 

the Newton -Raphs on iteration equation can be written: 


(45) 


3f 9f 

1u ' AU + ~W 


AV + 


9 ii 

3W 


AW = - f. 


(46) 


3£ 


2 

3U 


3f. 


AU + 


_ 2 _ 

3V 


3f. 


AV + 


2 

3W 


AW = - f- 


(47) 


3£, 3f_ 3£_ 

— — - • AU + • AV + 

8U av av 3W 


• AW = - f„ 


(48) 


where, for the i^ iteration. 


U. = U. 1 + AU 
i l-l 


(49) 


V. = V. + AV 
i l-l 


(50) 


W. = W. + AW (51) 

l i-l ^ J 

For an initial approximation U , V and W , the initial values of f, , f 0 

o o o 1 l 

and f^ from equations 42-44 and their associated partial derivatives 
(equation 45) can be calculated. From a matrix analysis of equations 
46-48, the values of AU, AV and AW can be determined and the values of 
U i s V i an< * W i ac *j uste< * according to equations 49-51. The above procedure 
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is repeated until a suitable convergence criteria is satisfied (i.e., 
until AU, AV and AW are less than some small number e) . 


Fluctuating Velocities, Turbulence Intensities and Stresses 
Equation 31 for hotwire sensor (1) is squared and time-averaged to 
give the following expression (neglecting small order terms such as 


4 2 2 

u , u v , etc.) : 


2 2 2 r 2 2 ,U.2 2 ,W. 2 _ r U. 

q L = a 5 v + [a 6 + a u + » 13 (y) + 2a 6 a u (y) 


„ ,W, „ ,UW. , 2 r 2 2 r W, 2 

+ 2a 6 a 13 M + 2a ll a l3 ( v 2 )] U + [a 7 a 12 M 


+ a 13 M + 2a 7 a l2 M + 2a 7 S 13 M + 2a 12 a l3 ^2^ W 


+ t2a 5 a 6 + 2a 5 a ll M + 2a 5 a 13 Ml UV 


+ [2a 5 a 7 + 2a 5 a 12 ( y ) + 2a g a 13 ( y ) ] vw 


+ [2a 6 a 7 + (2a 6 a 12 + 2a 7 a 13 ) v + (2a 6 a 13 + y 


+ (2a 11 a 12 + 2a 13 2 ) + 2a 11 a 13 ( v ) + 2a 12 a 13 ( y ) ] uw + 0(e) 


( 52 ) 


2 2 

Similar equations for q 2 and q 3 can be derived for hotwire sensors 

(2) and (3). 
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Likewise, equation 31 can be multiplied by the corresponding 
expression for and time -averaged to obtain: 


Il’z * *5 b s 1,2 * 


[a 6 b 6 * 3 ll b ll 


A 2 . ,W. 2 

¥ + a 13 b 13 ¥ 


rr w irur y 

+ ^ a 6 b ll + a ll b 6' ) V + (a 6 b 13 + a l3 b 6' ) V + (a ll b 13 + a 13 b ll^ U 


+ [a y b ? + a 12 b 12 ( y ) + a 1;J b 13 ( y ) + (a 7 b 12 + a 12 b ? ) y 


U UYi 2 

+ (a 7 b 13 + hsV V + (a l2 b 13 + ^ W + [a 5 b 6 + a 6 b 5 


+ Ca 5 b ll + a ll b 5^ 7 + (a 6 b 12 + a l2 b 5 ) V ] UV + [a S b 7 + a 7 b 5 


+ (a 5 b 12 + a 12 b S ) V + (a 5 b 13 + a l3 b 5 ) £ ™ + [a 6 b 7 + a 7 b 6 


+ C a 6 b 12 + a l2 b 6 + a 7 b 13 + a l3 b 7 > V + (a 6 b 13 + a 13 b 6 + a 7 b ll + a llV 


+ ( ' a il b 12 + a l2 b ll + 2a l3 b 13' ) y 2 + fa il b 13 + a l3 b l¥ y 2 


+ (a 12 b 13 + a 13 b 12 } UW + 


(53) 


where the constants b , n = 5 to 13 refer to the coefficients in the 

n' 

expression for q^ and are similar to the values of a^ in equation 28. 

Expressions similar to equation 53 can be derived for q"^q and 
q 2 q^. Thus, the set of equations 52 and 53 form a set of six non- 
linear equations in nine unknowns. As an approximate method of solution, 
the mean velocities U, V and W are first calculated from equations 39-41 


<|C 
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using the iteration method described in the previous section. These 

velocities are then used to calculate the coefficients of equations 52 

and 53. As a result, these equations are reduced to six linear equations 

in six unknowns which can then be solved simultaneously to give the 
~~2 ~~2 ~~2 

quantities u , v , w , uv, uw and vw. A computer program, coded in 
Fortran IV and given in Appendix C, has been written to use the hotwire 
data obtained from the experimental portion of this thesis to solve 
equations 39-41, 52 and 53 for the mean velocities, turbulence intensities 
and stresses respectively. The results are presented in the next 
chapter . 



II I II II III 
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EXPERIMENTAL RESULTS AND COMPARISON WITH THEORY 

Results of the experimental investigations are given in this 
chapter. In addition, comparison of these results with the theoretical 
results of the inviscid and viscid exact analysis programs are 
presented . 


Blade Static Pressure 

Experimental results are plotted in Figs. 12-16 for the five 
radial passage stations defined previously. It should be reiterated 
that the measurement stations do not correspond to constant radii, since 
the annulus passage is continuously varying. The measurement stations 
are illustrated in Fig. 3 and the pressure tap locations are 
specified in Tables 1 and 2. 

The inducer design characteristic of trailing edge loaded blades is 
apparent from the measured ^ distributions . measurements on the 
blade pressure surface remain positive across the entire chord length, 
with the gradient increasing continuously from hub to tip. The pressure 
surface ^ distribution decreases near the trailing edge, varying in the 
location at which the downswing begins from approximately 80% chord near 
the tip to greater than 90% chord near the hub. The blade suction 
surface measurements near the hub leading edge begin negative and 
become positive beyond 35% chord. At radial stations 3 thru 5, corres- 
ponding to mid-passage thru tip, the suction surface ^ distribution 
appears to begin with positive values, cross to negative values at 
approximately 20-30% chord and then return to positive values at 40-55% 
chord. The cross-over points increase in distance from the leading edge 
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Figure 13. Blade Static Pressure Distribution - Radial Station 2 
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at the radial stations nearer the tip. It is also noted that the 
difference between the pressure surface and suction surface at 
most chordwise locations increases continuously from hub to tip. 

It is apparent from Figs. 12-16 that there is a large discrepancy 
between the experimental and design curves. The design curves are based 
on two-dimensional theory and are derived from the mean streamline 
method of Wislicenus (Ref. 25) . The experimental results indicate that 
the three-dimensional inviscid effects are appreciable. 

The radial variation of the passage -averaged blade static pressure 
coefficients are compared in Figs. 17a-b for axial flow survey stations 
1 and 2 (shown in Fig. 3) . It can be seen that, for both locations, the 
discrepancy between design and experiment is greater near the tip, 
indicating the presence of velocity deficiencies resulting. from increased 
flow losses in this region. At station 1, the difference between the 
design and experimental curves does not appear to increase until 
approximately R = .8, whereas at station 2, the difference begins 

increasing at a radius much closer to the hub. This tends to indicate 
an increase in the extent of the loss region as the flow moves downstream 
through the inducer passage. 

Blade Limiting Streamline Angles 

The blade limiting streamline angle a is the limiting position (in 
degrees) of the flow streamline as the blade surface is approached 
(Fig. 2). The angle is measured from the two-dimensional or design 
flow direction and hence represents the extent of three-dimensionality 
in the flow. A similar parameter e w = tan a can also be defined which, 
in the peculiar geometry of the inducer, can approximate the ratio U/V 
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at the blade surface. The tangential variation of a with blade chord 
for each of the five radial measurement stations is shown for the 
pressure surface in Fig. 18 and for the suction surface in Fig. 19. 

The pressure surface distribution of a at the tip (radial station 
5) indicates negative values of a (and, thus, radially inward flow) from 
leading edge to mid -chord position. This tends to indicate the presence 
of the annulus wall boundary layer scraping effect which induces flow 
away from the tip. At all other radial stations, a increases 
continuously from leading edge to trailing edge. Near the hub trailing 
edge, a increases quite rapidly. The blade limiting streamline angles 
at both radial stations 1 and 2 appear to extrapolate beyond 90°, which 
is an indication of the existence of backflow in this region. This is 
presumably brought about by large radially outward flow that exists in 
the wake immediately downstream of the trailing edge. This has a 
tendency to decrease axial velocity near the hub and thus induce back- 
flows. At most axial locations, a decreases continuously from hub to 
tip. In several instances, this decrease appears linear. 

The suction surface a distribution remains relatively constant at 
all radial stations up to approximately 60% chord from the leading edge, 
when a more pronounced increase is noticed. At all stations except the 
tip, this increase extends to approximately 85% chord and then a 
decreases toward the trailing edge. This is possibly due to the blade 
blockage effect in this region. At the tip, a increases continuously 
and no decrease is noted. Again, as in the pressure surface distribution, 
a decreases continuously from hub to tip at practically all axial 
locations, and at some locations the variation appears linear. 



Limiting Streamline Angle a (degrees) 





73 


In most instances, the magnitudes of a on the suction surface are 
lower than at the corresponding position on the pressure surface. Fig. 
20 shows the radial variation of a at the trailing edge. The deviation 
between pressure and suction surface measurements decrease continuously 
with increasing radius from hub to tip. The magnitude of a, which is 
an indication of the extent of radial flows, is much higher than the 
values of a single blade reported in Ref. 14. This indicates that the 
radial velocity in the inducers are quite appreciable, especially near 
the blade surfaces. 


Mean Velocity Profiles 

The triaxial hotwire probe was used to measure the relative 
velocity profiles inside the inducer passage. As an indication of the 
effectiveness of this method in obtaining the relative velocity measure- 
ments, Fig, 21 compares the total relative velocity profile at station 1 
derived from the hotwire measurements with the results of Ref. 7 obtained 
from rotating pressure probe measurements. Good agreement is indicated 
at the two radii shown. 

The axial, radial and relative tangential velocity components 
described in this section are derived from the hotwire measurements and 
analysis described in the previous chapter. 

Measurements at Station 1 

Total Relative Velocity. - Fig. 22 shows the tangential variation 


of total relative velocity across the inducer passage at several 






Resultant Relative Velocity Q 


75 



0 10 20 


30 40 50 60 70 80 90 

Passage Width (% from suction surface) 


100 


Figure 21. Comparison of Rotating Hotwire and Rotating 
Pressure Probe Measurements 
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radii. A distinct velocity deficiency is noted at approximately 55% 
from the blade suction surface for all radial stations, but it is 
especially pronounced near the tip. This is the region of maximum loss. 
The radial flows inside the pressure and suction surface boundary 
layers, when encountered by the annulus wall, tend to roll toward mid- 
passage, interact, and produce strong vortices and radially inward 
flows. This interaction region is an area of considerable flow mixing, 
resulting in strong eddies and the associated energy dissipation. A 
concentration of high turbulence intensities in this region is confirmed 
by qualitative measurement of the A. C. fluctuating hotwire voltages. 

The radial variation of passage -averaged total relative velocity is 
plotted in Fig. 23. The difference between the design and experi- 
mental curves increases near the tip, further substantiating the 
existence of three-dimensional effects and flow loss in this region. 

The degradation in flow velocity near the tip also explains the 
behavior of the IjT variation in Fig. 17a. 

From the velocity profiles of Fig. 22, it is easy to discern the 
suction surface boundary layer at radii above R = .671. The suction 
surface boundary layer appears to grow in thickness as the tip is 
approached, increasing to approximately 25% of the passage width. This 
observation is consistent with the previous discussion about tip 
boundary layer interaction. No evidence of the pressure surface 
boundary layer can be detected in Fig. 22. This tends to indicate 
that the suction surface boundary layer is thicker than that of the 
pressure surface, although it should be remarked that no measurements 
were taken close to the blade surface. Since the blade element is not 
radial, the hotwire probe could not be located very close to the blade 
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surface while also assuring the correct orientation necessary for the 
data analysis described in the previous chapter. 

Relative Tangential Velocity .- Fig. 24 shows the passage 
variation of relative tangential velocity V at the several radii. The 
magnitude and shape of these curves are almost identical to the total 
relative velocity curves of Fig. 22, indicating the dominance of the 
tangential flow within the inducer passage. Comments concerning the 
total relative velocity are also applicable here. Fig. 25 is a plot 
of the radial variation of passage -averaged relative tangential velocity 
and indicates the region of large flow loss that exists near the tip. 

The absolute tangential velocity can be derived from 

C = R Q - V (54) 

u 

The high values of absolute tangential velocity near the tip indicates 
that the absolute stagnation pressure rise in this region is very large. 
This large absolute stagnation pressure rise is not due to flow turning 
but to the effects of complex viscous interactions. 

Axial Velocity . - Fig, 26 shows axial velocity W plotted versus 
percentage of passage width. The general trend for the tangential 
variation of axial velocity indicates an increase from suction surface 
to pressure surface. The radial variation of the axial velocity shows 
the largest values occurring near the hub, decreasing consistently 
towards the tip. This tends to indicate the effect of blade blockage 
on the axial velocity distribution. It is noted from Fig. 26 that 
negative values of W occur at the tip location R = .973. The existence 
of negative axial velocities at the extreme tip location indicates the 
presence of the annulus wall boundary layer scraping effect and was 
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similarly noted by the limiting streamline angle measurements at 
this location. 

It should be pointed out that the hotwire sensors used in the 
experimentation are not capable of distinguishing the direction of a 
velocity, only its magnitude. However, through the analysis of the 
previous chapter, hotwire equations are derived which assume an R-0-Z 
coordinate system and require appropriate direction cosines to the 
orientation coordinate system of the hotwire. The resulting system of 
equations are solved numerically. It is from the numerical solution 
of these equations that negative values for the velocity components can 
appear, indicating that the positive axis of that particular velocity 
component was actually 180° from that assumed in the measurement of the 
direction cosines. Thus it is possible, with the method derived in the 
previous chapter, to determine the magnitude and sense of the velocity 
vector measured by the hotwire sensors. 

The radial variation of passage -averaged axial velocity is shown in 
Fig. 27a. The rapid decrease in axial velocity near the tip is evident. 
The validity of the axial velocity distribution can be ascertained by 
applying the continuity equation to the experimental results. Using 

7 r t 

7T r t W x = / p2ir Wrdr (55) 

r h 

where is the uniform axial velocity upstream of the inducer, the 
results agree favorably with the design value. 

Radial Velocity .- Fig. 28 depicts the tangential variation of 
the radial velocity U. - Large values of U are found near the suction 
surface at radii close to the hub, indicating higher radially outward 



Inviscid analysis 

Viscid analysis 

Design 

• Experimental 





1 


J 













86 


flow in this region. Negative radial velocities are found at radii 
near the tip (R = .781 and greater), appearing at approximately 45% 
passage width. This is consistent with the previous discussions on 
boundary layer interaction and radially inward flow in this region. 

The rationale for obtaining negative velocity component values from the 
hotwire measurements was given in the previous section. Fig. 29a 
shows the radial variation of passage -averaged radial velocity. The 
values of U are quite large, indicating the appreciable three- 
dimensionality of the inducer flow. The radial velocities are of the 
same order of magnitude as the axial velocity. Fig. 29a indicates 
that the radial velocities are higher near the hub which confirms the 
conclusions of the blade, limiting streamline measurements. 

Measurements at Station 2 

Total Relative Velocity .- The tangential variation of total rela- 
tive velocity Q is shown in Fig. 30a for each of the measuring 
stations. Again, as in station 1, a region of distinct velocity 
deficiency is noted near the tip. The explaination for the velocity 
deficiency in this area has been attributed to the large flow losses 
encountered as the result of boundary layer interaction and extensive 
flow mixing discussed fully in an earlier section. In comparison with 
the results of station 1 (Fig. 22) , the position of the loss core 
appears to have shifted toward the suction surface to approximately 
40% passage width. A growth in the dimensions of the eddy 
inside the passage is evident as the flow proceeds from station 1 to 


station 2, 
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Figure 29. Radial Variation of Passage -Averaged 
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The radial variation of passage -averaged total relative velocity 
is plotted in Fig. 31. The values of Q R are considerably less than 
the design values, substantiating the effects of three-dimensionality 
on the inducer flow. The difference between the design curve and 
experimental results increases noticeably near the tip, further defining 
the region of large flow losses discussed previously. Comparison of 
Fig. 31 and Fig. 23 confirms the increase in size of the loss core 
within the blade passage as the flow proceeds downstream. This 
observation is also consistent with the behavior of the ip variations 
shown in Fig. 17. 

The extent of boundary layer growth on both the pressure and 
suction surface can be vaguely discerned in Fig. 30a, extending to 
approximately 20% passage width on the suction surface and to approxi- 
mately 10% passage width or less on the pressure surface. 

Relative Tangential Velocity. - Fig. 32a gives the variation of 
relative tangential velocity V across the passage width. Deviations 
from the total relative velocity profiles of Fig. 30a are slight, 
again indicating the dominance of tangential flow within the long, 
narrow inducer blade passages. In Fig. 33, the radial variation of 
passage -averaged relative tangential velocity is plotted. The region 
of large decrease in relative velocity can be easily discerned. 

Comments concerning total relative velocity in the previous section are 
also applicable here. This plot shows a significant departure from 
design values at all radii, the difference increasing rapidly as the 
tip is approached. These low relative velocities (or high absolute 
velocities) indicate an extremely large absolute stagnation pressure 
rise within the region and, thus, are an indication of the extent to 
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which complex viscid interactions are taking place. Comparison of 
Fig. 33 with Fig. 25 for station 1 clearly demonstrates the increase 
in size and severity of the loss core as the flow within the inducer 
passage proceeds downstream. 

Axial Velocity. - The tangential variation of axial velocity W is 
given in Fig. 30b for the various experimental radii. Overall magni- 
tudes are, of course, higher than those measured at station 1 (Fig. 26) 
due to the converging annulus. Again, as in station 1, the radial 
distribution of passage -averaged axial velocity (Fig* 27b) shows larger 
values occurring near the hub indicating the continuing presence of the 
blade blockage effect. It is interesting to note that the opposite 
trend was found in Ref. 20 at locations downstream of the trailing edge 
where no blade blockage effects should be present. This implies that 
significant changes occur in the axial velocity profile as the flow 
leaves the rotating inducer channel and proceeds downstream. These 
changes may be responsible, in part, for the backflow region previously 
reported in Ref. 10 near the hub trailing edge and confirmed by. the 
blade limiting streamline angle measurements of this report. Due to 
physical restraints, rotating hotwire measurements were not conducted 
close enough to the inducer hub to permit detection of backflows near 
the hub surface. The decrease in axial velocity near the tip, as shown 
in Fig. 27b, indicates the continuing presence of the annulus wall 
boundary layer scraping effect. The effect, however, is not as severe 
as at station 1 where negative axial velocities were measured (Fig. 27a) , 
Application of the continuity equation (equation 55) to the experimental 
results of Fig. 27b indicates good agreement with the design value. 
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An examination of Fig. 30b reveals a well-defined pressure 
surface boundary layer at R = .671 and R = .781 which extends approxi- 
mately 15% of the passage width. It appears that the boundary layer is 
thicker in this region than at the tip locations. Conversely, the 
suction surface boundary layer is well-defined near the tip at R = .973 
and R = .945, extending approximately 30% of the passage width and 
indicating that the suction surface boundary layer increases in thickness 
as the tip is approached. 

Radial Velocity .- The variation of radial velocity U across the 
passage is shown in Fig. 32b. The overall magnitudes appear larger than 
at station 1. Thus, the three-dimensional flow effects will be greater 
at station 2 and therefore accounts for the greater deviation of the 
flow from the two-dimensional design values which has been observed at 
this location (Fig. 31). Fig. 32b indicates negative radial velocities 
for the radii near the tip at approximately 25% from the suction surface. 
The radially inward flow at this location agrees with the previous 
discussions on boundary layer interaction and flow mixing which result 
in the velocity deficiencies and flow loss experienced in this region. 
Fig, 29b shows the radial variation of passage -averaged radial velocity. 
The significant radial flows are evident and this reflects the extent of 
three-dimensionality in the flow. As in station 1, the radial velocities 
are of the same order of magnitude as the axial velocities. Fig. 29b 
indicates that the radial velocities are higher near the hub, and 
comparison with Fig. 29a confirms the increase in radial velocities as 
the flow proceeds from station 1 to station 2, Both of these results 
confirm the observations of the blade limiting streamline angle measure- 
ments discussed in a previous section. 
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It is apparent from the discussions of radial flows at stations 1 
and 2 that radial velocities are significant within the inducer passages 
and, thus, any serious attempt to predict inducer flows must rely on a 
three-dimensional analysis. 


Turbulence Intensities and Stresses 
The distributions of tangential, axial and radial turbulence 
intensities at station 2, nondimensionalized with respect to local 
total relative velocity, are shown in isocontour form in Figs. 34, 35, 
and 36 respectively. The contours for all three intensity components are 
essentially the same, showing a "pocket" or "core" of high turbulence 
centered at approximately 40% passage width and R = .890. This coin- 

cides with the location of the maximum total relative velocity deficiency 
noted in Fig. 30a. The turbulence intensities are generally higher 
than those encountered in stationary passage. The peak intensities occur 
in the mixing region near the tip, where the two boundary layers merge 
and generate considerable flow mixing. The flow energy dissipated 
during this process is responsible for the velocity deficiencies 
encountered near the tip region in Fig. 30a. Another concentration of 
high turbulence is noted near the hub pressure surface and is an indi- 
cation of the proximity to the pressure surface boundary layer. The 
radial turbulence intensities ^ ^2 are generally higher than those in 
the axial and tangential directions, an indication of the violent 
radial motions occurring within the long narrow passages of the inducer. 
An isocontour plot of total turbulence energy at station 2, defined as 
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is given in Fig. 37. It essentially reflects the observations 
stated above, indicating high turbulence energy regions near the tip 
at approximately 40% passage width and near the hub pressure surface. 

It should be remarked that all the measurements reported here were 
taken away from the blade surfaces. The maximum turbulence intensities 
and stresses are likely to occur near the blade surfaces. High turbu- 
lence intensities measured away from the blade surfaces and reported 
here reflect the extent of turbulent mixing even near the mid-passage. 
The classical assumption that the viscous and turbulence effects are 
confined to very thin regions near the blade surfaces is evidently 
inapplicable to inducers and hence, a fully three-dimensional treatment 
is needed for the prediction of inducer flows. 

The distributions of turbulence velocity correlations uv", uw~ and 

2 

uw at station 2, nondimensionalized with respect to Q R , are given in 

isocontour form in Figs. 38, 39, and 40 respectively. These 
correlations are indicative of the stresses occurring within the 
inducer passage. Concentrations in stress intensities are similar to 
the turbulence intensity contours discussed previously. The maximum 
stresses occur in the mixing region near the tip at approximately 40% 
passage width. The radial stresses are by far the most dominant and 
emphasize the significant extent of three-dimensionality and complex 
viscous interaction occurring within the inducer channel, especially in 
the mixing region. As mentioned previously, the stresses shown in 
Figs. 38-40 represent values away from the blade surfaces. Stress 
values near the blade surface are likely to be high. The higher stress 
values noted near the hub pressure surface in Figs. 38-40 are an indi- 
cation of the proximity to the pressure surface boundary layer. 
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The distribution of turbulence stress/intensity ratios (uv/q ) , 

2 2 

(vw/q ) and (uw/q ) are shown in isocontour form in Figs. 41, 42, and 

2 

43 respectively. The magnitudes of uv/q vary from 0.01 in isolated 
points to 0.32 in the high loss mixing region. The contours generally 
follow those of uv" shown in Fig. 38, with ratios averaging approxi- 
mately 0.25 applicable in the regions of higher stress. The values of 
2 

vw/q range from 0.01 to 0.12, the higher magnitudes generally 

2 

occurring in the high stress areas. The distribution of uw/q indi- 
cates two regions where the ratios are high, corresponding to approxi- 
mately 40% and 65% passage width at R = .890. The magnitudes of 

2 

uw/q vary from 0.01 to 0.25, the higher values limited to the two 
regions defined above. 

Comparison of Experimental and Theoretical Results 
The Cooper -Bosch exact analysis program incorporating the modifi- 
cations mentioned in the chapter entitled "Theoretical Analysis" was 
run for the three -bladed Penn State inducer geometry. Both inviscid and 
viscid cases were considered. The inviscid program was run for 
approximately 75 relaxation cycles, resulting in a total RMS residual 
of 0.077. The viscid program was run for approximately 50 cycles and 
produced a total RMS residual of 0.200. In this section, the results 
of the inviscid and viscid analyses will be discussed and compared with 
the experimental results described earlier in this chapter. 
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Blade Static Pressure 


The inducer blade static pressure distributions derived from the 
exact inviscid analysis are plotted in Figs. 12-16 at the locations 
corresponding to the experimental ij'g measuring stations. As can be 
seen, the results agree remarkably well at all radial stations and 
especially for those near the hub (radial stations 1 and 2). Close to 
the hub, where three-dimensional viscid effects are not prevalent, the 
difference between the experimental and theoretical results should not 
be large. Near the tip region, the experimental pressure surface \p s 
distributions agree closely with the numerical analysis while the 
suction surface \p s distributions show the most discrepancy. The radial 
variation of ip s plotted in Fig. 17 shows the larger deviation between 
the experimental and theoretical results near the tip which is expected 
from previous discussions and reflects the region where the secondary 
flow effects are concentrated. The chordwise gradients of pressure and 
suction surface i^ s in Figs. 12-16 appear similar for both experimental 
and theoretical results . 

The distributions obtained from the viscid analysis program are 
also plotted in Figs. 12-16. The general shape of the viscid analysis 
4> s distribution closely resembles that of the inviscid analysis. 

Greatest variation between viscid and inviscid analyses are shown at 
the tip (radial station 5), where the viscous effects are dominant. 

The radial variations of \ in Fig. 17 verify this observation. The 
deviation between inviscid and viscid near the tip appears greater 
at station 2 than at station 1. 
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Mean Velocities at Station 1 

Total Relative Velocity .- Hie inviscid analysis results for total 
relative velocity Q R across the inducer passage are plotted in Fig. 

44 for the radii corresponding to those used in the experimental investi- 
gation. The similarity of the analytical solution with experiment can 
be seen by comparing the theoretical results with Fig, 22. As in the 
experimental results, the velocities increase from hub to tip. Inviscid 
turning effects resulting in higher velocity near the suction surface 
are also evident. Overall magnitudes of the velocities are lower than 
those found from experimentation. Thus, although the static pressure 
distributions have been found to be comparable, the velocities 
predicted are lower than the measured values and seems to indicate loss 
in kinetic energy. The radial variation of passage -averaged total 
relative velocity Q found from the inviscid analysis is compared in 
Fig. 23. The difference between the inviscid results and design 
values of appears to increase from hub to tip. This tends to indi- 
cate that three-dimensional inviscid effects are larger near the tip, 
a result which has been shown experimentally. 

The total relative velocity distribution obtained from the viscid 
analysis is shown in Fig. 45. The imposition of the boundary condi- 
tion which defines the relative velocity on the inducer blade surface 
as zero enables the viscid program to provide a crude approximation for 
the pressure and suction surface boundary layers. The magnitudes of 
Q r are similar to those found from the inviscid analysis. A slight 
velocity deficiency is noted near the tip at approximately 50% 
passage width. This agrees with the experimental results of Fig. 22 
and indicates an area of high viscous loss. A closer comparison of 
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Figure 44. Inviscid Analysis Results for Total Relative Velocity 

Distribution at Station 1 
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results can be made by examining Fig. 46 for R = .973 and Fig. 47 for 
R = .548. The viscid results seem to agree with experiment qualita- 
tively. TTie radial variation of Q|^ in Fig. 23 shows that, near the tip, 
the viscid analysis velocities are lower than those derived from the 
inviscid analysis. 

Relative Tangential Velocity. - The relative tangential velocity 
distributions obtained from the inviscid and viscid analyses are not 
included here, since the magnitude and shape of these curves closely 
approximate the total relative velocity distributions of Figs. 44 and 
45 and illustrates the dominance of the tangential component inside the 
inducer channel . 

Axial Velocity .- The axial velocities predicted from viscid and 
inviscid analysis at R = .973 and R = .548 are shown compared with 
experimental results in Fig. 48. The predictions were found to be 
very similar at all radii. The predictions are good at R = .548 and 
poor at R = .973, thus indicating the dominance of viscous effects near 
the tip. No comments can be made with regard to the accuracy of viscid 
results, since the measurements close to the proximity of the wall are 
not available. The passage -averaged velocities W are compared with 
experimental results in Fig. 27. The predictions are good up to R = .9, 
the discrepancy increases considerably beyond this radius. 

Radial Velocity .- The radial velocity predicted from the theoreti- 
cal analyses is found to be very small at this station. The radial 
variation of passage -averaged radial velocity U, plotted in Fig. 29, 
indicates that the theoretical predictions are significantly lower than 
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Figure 46. Comparison of Total Relative Velocity 
Distributions at Station 1, R = .973 
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Figure 47. Comparison of Total Relative Velocity 
Distributions at Station 1, R = .548 
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Figure 48. Comparison of Axial Velocity Distributions 
at Station 1, R = .973 and R = .548 
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the experimental results and emphasizing the substantial three- 
dimensionality which exists in the real inducer flow. 

Mean Velocities at Station 2 

Total Relative Velocity. - Inviscid analysis results for the total 
relative velocity distribution are plotted in Fig. 49. Magnitudes of 
Q near the tip are comparable to those found experimentally (Fig. 30a), 
whereas near the hub the inviscid velocities are significantly lower 
than those shown for the experimental results of Fig. 30a. The latter 
observation is consistent with the results for the Q R distribution at 
station 1. The radial variation of passage -averaged total relative 
velocity predicted from the inviscid analysis is shown in Fig. 31. The 
difference between the inviscid results and design values of Q appears 
to increase from hub to tip, indicating the larger three-dimensional 
inviscid effects existing near the tip and confirmed by experimentation. 
Comparison of inviscid results from Fig. 31 with those of Fig. 23 for 
station 1 indicates that the deviation between the inviscid and design 
velocities has increased from station 1 to station 2. This supports 
the contention that the three-dimensional inviscid effects increase in 
severity as the flow proceeds downstream inside the inducer channel. 

The agreement between the measured and predicted Q R is reasonably good 
(Fig. 31). 

The viscid analysis prediction for the total relative velocity 
distribution is shown in Fig. 50. It is a striking departure from the 
inviscid analysis distribution (Fig. 49), especially near the tip 
where the viscous effects are known to be appreciable. The viscid 
analysis also provides crude approximations for the suction and pressure 
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R = .973 

R = .945 

R = .890 



Figure 49. Inviscid Analysis Results for Total Relative 
Velocity Distribution at Station 2 
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surface boundary layers. A large velocity deficiency near the tip is 
predicted at approximately 50% passage width and agrees favorably with 
the experimental Q R profiles plotted in Fig. 30a. It is apparent 
that the velocity deficiency noted near the tip at station 1 (Fig. 

45) has grown considerably as the flow proceeded downstream to 
station 2, indicating an increase in size and intensity of the viscous 
loss region and substantiating the experimental results discussed 
previously. The blade boundary layer development predicted by the 
viscid analysis can be seen in Fig. 50. The suction surface boundary 
layer appears thicker than the pressure surface boundary layer at all 
radii, increasing in thickness from hub to tip. This observation is 
consistent with experimental results and conforms with previous 
discussions on boundary layer interaction and flow mixing near the tip. 
The radial variation of Q R derived from the viscid analysis is plotted 
in Fig. 31. The deviation between viscid and inviscid velocities near 
the tip can be attributed to the large viscous losses which are known 
to exist in this region. Comparisons of the total relative velocity 
distributions at R = .973 and R = .548 are given in Fig. 51 and 
Fig. 52 respectively. It appears that the viscid analysis distri- 
bution provides better approximations to the experimentally derived 
velocity profiles. 

Relative Tangential Velocity .- Relative tangential velocity distri- 
butions predicted by the inviscid and viscid analyses closely resemble 
the total relative velocity distributions indicated in Fig. 49 and 
50 respectively. Comments in the previous section relating to total 
relative velocity are also applicable here. The tangential velocity 
continues to dominate the flow within the rotating inducer passage. 
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Figure 51. Comparison of Total Relative Velocity Distributions 

at Station 2, R = .973 
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Figure 52. Comparison of Total Relative Velocity 
Distributions at Station 2, R = .671 
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although the influence of radial and axial velocities is greater than 
that at station 1, especially near the hub. The three-dimensionality 
of the flow is appreciable. The radial variation of passage -averaged 
relative tangential velocity derived from the inviscid analysis is 
shown in Fig. 33. The influence of three-dimensional inviscid 
effects on the flow, reflected in the difference between inviscid 
results and design, appears to increase in magnitude near the tip. The 
high absolute tangential velocities predicted from the inviscid 
analysis are indicative of the high stagnation pressure rise through 
the inducer passage. The lower values of viscid analysis V, when 
compared to those predicted from the inviscid analysis (Fig. 33), can 
be attributed to the viscous losses and secondary flows which prevail 
at this location. The higher values of absolute tangential velocity 
predicted from the viscid analysis are consistent with the increased 
effects of complex viscous interactions near the tip. 

Axial Velocity .- Fig. 53 shows the axial velocity distribution 
predicted from the inviscid analysis. There is a definite decrease in 
axial velocity from hub to tip, which appears to indicate the presence 
of the blade blockage effect within the flow passage. This observation 
is consistent with the experimental results plotted in Fig. 30b. The 
radial variation of W predicted from the inviscid analysis (Fig. 27) 
agrees almost exactly with the experimental distributions . Continuity 
has been satisfied within the exact analysis program. 

The viscid analysis results for the axial velocity distributions 
are shown in Fig. 54 and indicate the approximate profiles for the 
pressure and suction surface boundary layers. The axial velocity 
profile decreases in magnitude from hub to tip and tends to confirm 
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Figure 53. Inviscid Analysis Results for Axial Velocity 
Distribution at Station 2 
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Figure 54. Viscid Analysis Results for Axial Velocity Distribution 

at Station 2 
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the presence of the blade blockage effect at this location. The radial 
variation of viscid analysis W in Fig. 27b conforms almost precisely 
with the experimental results. Comparisons of the axial velocity 
distributions at R = .973 and R = .671 are given in Fig. 55. The 
axial velocity profiles predicted from the viscid analysis appear to 
more closely approximate the experimental distributions . 

Radial Velocity . - The inviscid analysis results for the radial 
velocity distribution at station 2 are given in Fig. 56. The tip 
region appears to exhibit the lowest radial velocities, which is 
consistent with the experimental results plotted in Fig. 32b. The 
blade limiting streamline angle measurements also substantiate this 
observation. The magnitudes of the inviscid analysis radial velocities 
are significantly lower than the corresponding values of experimental 
radial velocity, indicating the considerable three-dimensionality of 
the real flow existing within the inducer blade passages. The radial 
velocity plotted in Fig. 56 appears to decrease across the inducer 
passage from suction surface to pressure surface, a condition which is 
found to exist experimentally (Fig. 32b). A region of radially 
inward flow near the tip at approximately 50% passage width is noted 
from the inviscid analysis radial velocity profiles. The existence 
of radially inward flow in this area has been found experimentally and 
is consistent with previous discussions on flow mixing effects in this 
region. 

The viscid analysis results for the radial velocity distribution is 
shown in Fig. 57. Little difference is noted from the inviscid 
distribution except at the pressure and suction surfaces where the 
velocities are fixed at zero. The region of radially inward flow near 
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Figure 56. Inviscid Analysis Results for Radial Velocity 
Distribution at Station 2 






129 


the tip is similar to that found with the inviscid analysis. The radial 
variation of passage -averaged radial velocity derived from the viscid 
analysis (rig. 29) reflects the significantly higher radial velocities 
found experimentally. 
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DISCUSSIONS AND CONCLUSIONS 

A brief summary of conclusions reached by the analytical and experi- 
mental investigation reported in this thesis are as follows: 

1) Preliminary investigations into the use of the quasi-three- 
dimensional Douglas -Neumann method indicate that it appears to provide a 
better initialization of velocity and pressure parameters needed for the 
Cooper-Bosch exact solution. 

2) Modifications to the Cooper-Bosch program to automatically 
unload the trailing edge station and incorporate dominant viscid effects 
have been made in subroutines "Main", "Dloss u , and "Resid". Preliminary 
testing of these modifications indicate that the inviscid analysis has 
been improved and a satisfactory viscous capability has been provided. 

3) The viscid analysis is, at best, approximate due to the various 
assumptions and simplifications made. In particular, the viscid boundary 
conditions imposed on the solution are rather drastic, since the grid 
geometry spacing used in the exact analysis is relatively large. More 
tangential grid stations would be needed, especially close to the blade 
surface, to better define the shape of the blade boundary layer. 

4) The ammonia trace technique provides a satisfactory method 

for determining blade limiting streamline angles within the rotating 
inducer blade passages. The blade limiting streamline angle measure- 
ments provide several observations which have either been noted in 
previous investigations or have been found from other experimental 
results contained in this thesis. These include: an increase in a 

from the leading edge to the trailing edge indicating the existence of 
higher radial velocities as the flow proceeds downstream within the 
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inducer channel; higher values of a near the hub indicating higher 
radial velocities in this region; negative values of a at the pressure 
surface tip (up to 45% chord) indicating radially inward flow due to 
the presence of the annulus wall boundary layer scraping effect; values 
of a greater than 90° near the hub trailing edge indicating the 
existence of a backflow region in this area; higher values of a through- 
out the inducer channel in comparison with the results of a single blade 
reported in Ref. 14 and indicating appreciable radial velocities 
existing within the inducer passage, especially near the blade surface. 
In addition, the blade limiting streamline measurements will provide 
valuable information on boundary layer and Reynolds stress character- 
istics for future investigations. 

5) The experimental blade static pressure distributions confirm 
the trailing edge loading characteristic inherent in the inducer blade 
design. The magnitudes of are considerably higher than design 
values, indicating the significant effect of three-dimensionality in 
the inducer flow. The static pressure distributions display small 
negative ^ values near the leading edge of the suction surface. The 
agreement between theory and experiment is good, especially near the 
hub where flow mixing and viscous loss effects are minimal. 

6) The rotating triaxial hotwire probe utilized in this study has 
yielded satisfactory velocity profiles and turbulence quantities. 
Comparison of velocities derived from this method show good agreement 
with those found from rotating pressure probe measurements in Ref. 7. 

It can be concluded that hotwire anemometry can be an extremely useful 
tool in the study of the relative flow parameters in a rotating 


environment. 
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7) Total relative velocity measurements indicate a substantial 
velocity deficiency near the tip at mid-passage which expands 
significantly as the flow proceeds downstream toward the inducer 
trailing edge. This indicates the considerable extent of flow mixing 
due to boundary layer interaction, radially inward and outward flow, 
annulus wall effect, etc. which is prevalent at this location inside 
the long narrow passages of the inducer. The position of this "loss 
core" appears to drift closer to the blade suction surface as the flow 
proceeds toward the inducer exit. An increased difference between the 
experimental results and the two-dimensional design curves near the tip 
for the radial distribution of passage -averaged total relative velocity 
further indicates the regions where significant three-dimensional 
effects and flow losses exist. The velocity distributions determined 
from the theoretical analysis are similar to those derived from experi- 
mentation. The presence of the velocity deficiencies near the tip has 
been predicted by the viscid analysis program. 

8) The high values of absolute tangential velocity which are 
found to exist near the inducer tip indicate a region of large 
absolute stagnation pressure rise caused by the effects of complex 
viscous interactions. The size and severity of this region increases 
significantly as the flow proceeds downstream inside the inducer blade 
channel. This observation is also predicted from the three- 
dimensional theoretical analysis. 

9) The higher values of axial velocity near the hub indicate 
the significant effect of blade blockage within the inducer flow 
passage. A slight backflow was found to exist at the extreme tip 
location of flow station 1 and can be attributed to the annulus wall 
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boundary layer scraping effect identified in the measurements . 

The agreement between the experimental results and theoretical analysis 
is extremely good. The viscid analysis results appear to more closely 
approximate the experimental axial velocity distributions than the 
inviscid program. 

10) The measured magnitudes of radial velocity are found to be on 
the same order as those of axial velocity within the inducer passage. 
This is similar to the results of Ref. 20 whicfr were made slightly 
downstream of the inducer trailing edge. The large values of radial 
velocity confirm the highly three-dimensional characteristic of inducer 
flow and emphasize the necessity of a suitable three-dimensional theory 
for accurate flow analysis. The radial velocities are generally found 
to decrease in magnitude from hub to tip and increase in magnitude from 
leading edge to trailing edge. These observations are consistent with 
the measurements mentioned previously. A region of radially inward 
flow is found to exist near the mid-passages of the tip region and 
correspond to the locations of the velocity deficiencies noted in the 

distributions. These measurements support the explanation for the 
existence of the large flow losses in this region due to extensive flow 
mixing and complex viscous interations. The magnitudes of radial 
velocity predicted from the theoretical analysis are significantly 
lower than the experimental results. The radially inward flow found 
experimentally near the tip mid-passage has been predicted with the 
three-dimensional theoretical analysis. 

11) Turbulence levels within the blade passage, indicated from the 
experimental results of this thesis, are generally high near the tip 
regions. A growing core of high turbulence is evident near the tip 
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mid-passage corresponding to the region of velocity deficiency 
mentioned in item 7) above. High turbulence levels are also recorded 
near the hub pressure surface and indicate the influence of the pres- 
sure surface boundary layer. The radial component of turbulence 
intensities appears to have the largest magnitudes, reaching calculated 
values of up to 24% in the mixing region. The total turbulence energy 
contours are similar to those of the turbulence intensities. 

12) The locations of high turbulence stresses are concentrated in 
the high turbulence intensity areas of the mixing region and near the 
hub pressure surface. Values of uv appear to be higher in the mixing 
region than the corresponding values of vw and uw. The high stress 
regions are indications of areas subjected to complex viscous inter- 
actions . 

13) Experimental velocity and turbulence results of Ref. 20 
show nearly uniform distributions slightly downstream of the inducer 
trailing edge. Thus it appears that considerable wake diffusion occurs 
immediately after the trailing edge to decay the turbulence core and 
blade blockage effects reported in this thesis. 

The hotwire analysis developed in this investigation was an 
initial attempt to determine the three components of velocity, 
turbulence intensity and Reynolds stresses in a rotating reference 
frame. Certain improvements and refinements can be made to the 
experimental techniques to reduce the errors encountered in the experi- 
mental measurement. In particular, the following are recommended: 

a) The use of linearizing circuits in conjunction with the hotwire 
anemometer for measurement of flow parameters when turbulence levels are 
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excessively high (greater than 20%) . The linearizing assumptions 
used for the hotwire equation derivation would then be valid. 

b) The use of a multiplying circuit instead of the sum -and - 
difference circuit to obtain the direct measurement of e^e 2 terms for 
the hotwire equations. This would minimize the experimental errors 
encountered in the hotwire signal handling. 

c) The necessity of having identical resistances for all hotwire 
sensors of the probe so the same calibration curve applies to all. This 
is extremely important in fluctuation velocity measurements where 
arithmetic operation of hotwire signals is involved. 

d) The necessity of precise measurement of hotwire angles with 
respect to the R-0-Z coordinate system. The velocities and turbulence 
quantities calculated by the analysis described in the chapter "Deriva- 
tion of the Hotwire Equations" are sensitive to the direction cosine 
coefficients used in the equations. 

e) The use of a precise traversing mechanism to survey the flow 
field and assure uniform orientation of the hotwire probe in the R-6-Z 
coordinate system. 

This thesis has attempted to present a detailed description of 
analytical and experimental investigations on flow through a three - 
bladed axial flow inducer. Although the investigations are performed 
on an axial flow inducer of a specific configuration, certain methods 
and techniques which have been developed are applicable to all types of 
turbomachinery. In particular, the equations and method of solution 
used in the exact analysis program are completely general and are not 
restricted solely to the solution of inducer fluid flow. Similarly, the 
experimental techniques used in conjunction with the rotating three- 
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sensor hotwire anemometer can be applied to most turbomachinery 
applications where relative velocity measurements are desired. 
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APPENDIX A 

Flow Chart Diagrams of Exact Analysis Program 
Incorporating Viscid and Inviscid Modifications 
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APPENDIX B 


Fortran Listing of Exact Analysis Program 
Incorporating Viscid and Inviscid Modifications 


c 

C THREE DIMENSIONAL FLOW ANALYSIS IN PENN STATE INDUCER JULY 1973 
C 

C ORIGINAL PROGRAM DEVELOPED BY P. COOPER AND H. BOSCH OF 
C TRW ACCESSORIES DIVISION, CLEVELAND, OHIO 

C 

C MODIFICATIONS MADE BY C. GORTON AT PENN STATE DURING SUMMER 1973 
C 

IMPLICIT LOGICAL*l ($) 

COMMON U(7,7,39),V(7,7,39),W(7,7,39),P(7,7,39),D(7,7,39),R(7,7,39) 
l,T(7,7 t 39),Z(7,7,39),FR,FT,FZ,DMAX(U),DELX(M,NSEQ(U),RES(8,8,U0), 
2AR(7,7,39),AT(7,7,39).AZ(7,7,39),BR(7,7,39),BT(7,7,39),BZ(7,7,39), 
3CR(7,7 ,39) ,CT(7,7 ,39) ,CZ(7 »7 *39) ,XH,RSTAR, 

l*TT,PSAT,DLIQ,REV t A,M,NVAR,NTR,IRSTAR,DX,I , j,k,ii, jj,kk,imax,jmax,k 
5MAX , Vise , NBD ,KLE ,KTE ,WRE (7 ,7 ,39 ) 

COMMON /$$LOSS/ TAUR(7,7,26) ,TAUT(7 ,7,26) ,TAUZ(7,7,26) ,TAU(7,7,26) 
COMMON / $ $ CFRE/ SLOPEP ( 25 ) , SLOPES ( 2 5 ) ,BP( 25 ) ,BS ( 2 5 ) 

COMMON /$$STRT/ KSTART 
DIMENSION CFS(2),CFP(2) 

DIMENSION THETA(3 ) ,REREF1 (3) ,REREF2(3) ,CFREFP(2,3) ,CFREFS( 2,3 ) 
DIMENSION DEV(7),DEL(7) 

READ INITIAL DATA FROM TAPE 

READ ( 91 ) IMAX, JMAX ,KMAX ,KLE , KTE 

READ(91)(((R(I,J,K) ,1=1, IMAX ),J=1, JMAX) ,K=1, KMAX ) 

READ(91 ) ( ( (T(I ,J,K) ,1=1 ,IMAX) ,J=1 , JMAX) ,K=1 .KMAX) 

READ ( 91 ) ( ( ( Z ( I , J , K) , 1=1 , IMAX ) , J=1 , JMAX ) ,K=1 , KMAX ) 

READ ( 91 ) KOUNT ,NUM,M, (NSEQ(L) ,]>1 ,1+ ) 

READ ( 91 ) ( Df 1AX ( L ) , L= 1 , U ) , A , CR I T , E 
READ ( 91 ) REV ,DLIQ ,PS AT , TT , VISC 

READ(91 )( ( (U(l,J,K) ,1=1, IMAX) ,J=1,JMAX) ,K= 1 ,KMAX) 

READ ( 9 1 ) ( ( ( V ( I , J , K ) , I = 1 , IMAX ) , J = 1 , JI 4AX ) , K=1 , KMA X ) 

READ (91 ) ( ( ( W( I ,J ,K) ,1=1, IMAX), J=l, JMAX), K=l, KMAX) 

READ(91)( ((P(I,J,K), 1=1, IMAX), J=l, JMAX), K=l, KMAX) 

REWIND 91 


JPRES=1 

JSUC=JMAX 

NBD=1 

NTR=0 

ANR=1* * IMAX* JMAX* KMAX 
UPDATE TAPE DATA 

READ( 5,1001+ )KOUNT,NUM,M, (NSEQ(L) ,L=1,1+) 

READ( 5 , 1002 ) ( DMAX ( L ) , L=1 , 1+ ) , A , CRIT ,E 
READ ( 5 ,1002 )REV ,DLIQ , PS AT ,TT , VISC 
IF(VISC)32,31,32 

IF VISCOSITY IS CONSIDERED, SET U,V,W EQUAL TO ZERO ON BLADE 


32 DO 17 K=KLE,KTE 
DO 17 1=1, IMAX 
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DO 17 J=1 , JMAX - 

IF( ( J-l )* ( JMAX- J ) ) 31 ,18 ,17 

18 U(I,J,K)=0.0 
V(I,J,K)=0.0 
W(I,J,K)=0.0 

17 CONTINUE 

CALCULATE INITIAL OUTLET AIR ANGLES 

31 DO 121 1=1 ,IMAX 
IE( VISC )19 , 20 ,19 

20 DEL(l)=ATAN(V(l,JFRES,KTE)/W(l,JPRES,KTE)) 

GO TO 21 

19 DEL(l) =ATAN( V ( I , JPRES , KTE+1 ) /W( I , JPRES ,KTE+1 ) ) 

21 DEL(l)=ABS(DEL(l)) 

DDEL=DEL( I )*l80 . /3 .14159 

121 CONTINUE 

CONVERT COORDINATE SYSTEMS 

KSTART=1 
CALL JACOB 
PRINT 1000 
NMAX=KOUNT+NUM 

PRINT 1007,IMAX,JMAX,KMAX,M,A,REV,DLIQ,PSAT,TT,VISC 
PRINT 1006,(NSEQ(L),L=1,4) 

READ IN CF VS. RE CURVES FOR VARIOUS REFERENCE THETA LOCATIONS 

IF ( VISC ) 669 , 67 0 , 669 
669 READ(5,700)NCURVE 
DO 710 1=1 ,NCURVE 

READ ( 5,720 )THETA( I) ,REREFl(l) ,CFREFP(l,l) ,CFREFS(l,l) ,REREF2(l) ,CF 
1REFP(2,I) ,CFREFS(2,I) 

THETA(l)=THETA(l)*3.l4l59/l80. 

710 CONTINUE 

DEVELOP EQUATIONS OF CF VS. RE CURVES FOR ALL BLADE K VALUES 

DO 750 K=KLE,KTE 
I=IMAX 

DO 755 Jl-1,2 
IF(J1.EQ.2)G0 TO 756 

Jl=l PRESSURE SURFACE Jl=2 SUCTION SURFACE 

J=1 

GO TO 757 

756 J=JMAX 

757 THET=ABS(T(I , J,K)-T(I , J,KLE) ) 

SEARCH THETA REFERENCE VALUES 
DO 760 KK=1,NCURVE 
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II=KK 

IF(THET.GT. THETA ( II) )G0 TO 759 
IP(II.EQ.1)G0 TO 762 
GO TO 761 
762 11 = 11+1 

761 ratio=(thet-theta(ii-i))/(theta(ii)-theta(ii-i)) 

GO TO 765 

759 IF(lI.EQ.NCURVE)GO TO 76l 

760 CONTINUE 

765 DO 7*»0 JJ=1,2 

GO TO (780,785) ,J1 

780 CFP(JJ)=CFREFP(JJ,II-1)+RATI0«(CFREFP(JJ,II)-CFREFP( JJ,II-l) ) 
GO TO 7^0 

785 CFS( JJ )=CFREFS( JJ ,11-1 )+RATI0*(CFREFS( JJ,Il)-CFREFS( JJ,II-1 ) ) 
7^0 CONTINUE 
755 CONTINUE 
JJ=1 

DEN0M= AL0G1 0 ( REREF2 ( J J ) ) -AL0G10 ( REREF1 ( J J ) ) 

SLOPEP (K )= ( AL0G10 ( CFP( JJ) ) -ALOGIO ( CFP( JJ+1 ) ) ) /DENOM 
SLOPES (K)=( ALOGIO (CFS(JJ)) -ALOGIO ( CFS ( JJ+1 ) ) ) /DENOM 
BP( K) =CFP( JJ ) *REREF1 ( JJ )**SLOPEP(K) 

BS ( K ) =CFS ( J J ) «REREF1 ( J J ) *»SLOPES ( K ) 

750 CONTINUE 

CALCULATE POINT DENSITIES 

INITIALIZE VALUES OF TAUR, TAUT, AND TAUZ IF NECESSARY 

670 DO 150 L=1,U 
150 DELX( L ) =DMAX ( L ) 

DO 501 K=1,KMAX 
DO 501 J=1,JMAX 
DO 501 1=1 , IMAX 
IF (TT) 502,503,502 

502 CALL STATE 
GO TO 120 

503 D(I,J,K)=DLIQ 

120 IF(VISC)385,501,385 
385 IF((K.LT.KLE).OR„(K.GT.KTE))GO to 501 
CALL DLOSS 
501 CONTINUE 

USE FRICTION VALUES TO FIND INITIAL FR , FT , FZ , AND RESIDUALS 

DO lUO K=1,KMAX 
DO lUO J=1,JMAX 
DO lUO 1=1, IMAX 
IF (D(I,J,K)) 500,500,500 

500 CALL RESID 
140 CONTINUE 
PRINT 1012 

CALCULATE TOTAL ROOT-MEAN-SQUARE RESIDUAL 
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467 RT=0.0 

RESMAX=0 . 0 
RESMIN=RES( 1,1,1) 

DO 368 K=1,KMAX 
DO 368 J=1 , JMAX 
DO 368 1=1 , IMAX 

IF(RES(l,J,K)-RESMAX)403,403,402 

402 RESMAX=RES(l,J,K) 

GO TO 368 

403 IF(RES(I,J,K)-RESMIN)404,368,368 

404 RESMIN=RES(l,J,K) 

368 RT=RT+RES(l,J,K) 

RMS=SQRT ( RT/ ANR ) 

RESMAX=SQRT ( RESMAX/4 . ) 

RESMIN=SQRT(RESMIN/4. ) 

PRINT 1013 .FOUNT, NTR, RMS, RESMAX.RESMIN, (DMAX(L) ,L=1,4 ) 

DO 151 L=1 ,4 
151 DMAX(L)=0 „ 0 

IF ( RESMAX-CRIT*E ) 369,369,370 

370 K0UNT=K0UNT+1 
NTR=0 

IF (KOUNT-NMAX) 360,360,371 
START SUCCESSIVE VARIATIONS CYCLE 

CHANGE OUTLET AIR ANGLE AND DOWNSTREAM COORDINATES IF NECESSARY 

360 KBEGIN=2 
KEND=KTE 
NTR=0 

$DELTA=. FALSE. 

381 DO 460 KK=KBEGIN,KEND 
DO 460 JJ=1 , JMAX 
DO 460 I 1=1, IMAX 
IRSTAR=1 
CALL STAR 
DO 460 L=l,4 
nvar=nseq(l) 

GO TO (111,222,333,450) ,NVAR 
111 IF (KK-2) 460,460,611 

611 IF ( (H-l)*(lMAX-II) ) 460,460,450 

222 IF (KK-2) 460,460,481 

481 IF (KK-KMAX) 450,482,460 

482 IF ( ( JJ-1 )*( JMAX-JJ ) ) 460,460,450 

333 IF (KK-2) 460,485,486 

485 IF ( ( JJ-l )*(JMAX-JJ ) ) 460,460,471 

471 IF ( (II-1)*(IMAX-II)) 460,472,450 

472 IF (AZ(II,JJ,KK)) 460,450,460 

486 IF (KK-KMAX) 487,460,460 

487 IF( (JJ-1)*(JMAX-JJ))460, 460,450 
450 CALL ADJ 

460 CONTINUE 

if(kend-kmax) 46 i , 152,152 

461 K2=KTE+1 
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DO 510 I= 1 ,IMAX 

* 

CHANGE EXIT AIR ANGLE IF GREATER THAN TWO PERCENT OF PREVIOUS 
EXIT AIR ANGLE 

PCNT2= . 02*DEL( I ) 

IF(VISC)25, 26,25 

26 DELNEW= ABS ( ATAN( V ( I , JPRES , KTE ) / W ( I , JPRES , KTE ) ) ) 

GO TO 27 

25 DELNEW=ABS ( ATAN ( V ( I , JPRES , KTE+1 ) / W ( I , JPRES ,KTE+1 ) ) ) 

27 ABDEL= ABS ( DELNEW-DEL ( I ) ) 

IF ( ABDEL-PCNT2 ) 51 0 , 54 0 , 54 0 

540 DEL(I)=DELNEW 

DDEL=DEL(l)*l80./3.l4l59 
$DELTA= . TRUE . 

WRITE(6,55)I,DDEL 
DO 515 K=K2 , KMAX 
DO 515 J=1 .JMAX 

T(l,J,K)=T(l,J,KTE)+(Z(l,J,KTE)-Z(l,J,K))*TAN(DEL(l))/R(l,J,K) 
515 CONTINUE 
510 CONTINUE 

RE-CALCULATE JACOBIAN COEFFICIENTS FOR COORDINATE TRANSFORMATION 

IF(.NOT.$DELTA)GO TO 521 

KSTART=KTE 

CALL JACOB 

RECALCULATE RESIDUALS 
NBD=1 

KKTE=KTE-1 
DO 520 K=KKTE,KMAX 
DO 520 J=1,JMAX 
DO 520 1=1 , IMAX 
CALL RESID 

520 CONTINUE 

521 KBEGIN=KTE 
KEND=KMAX 
GO TO 381 

152 DO 153 L=l,4 

IF (DMAX(L) ) 561,562,561 

561 DELX(L)=DMAX(L) 

GO TO 153 

562 DELX(L)=DELX(L)*A 

153 CONTINUE 
GO TO 467 

*** OUTPUT ROUTINE *** 

98 WRITE(92 )D1AX, JMAX, KMAX, KLE, KTE 

WRITE(92) ( ( (R(I ,J,K) ,1=1, IMAX) ,J=1,JMAX) ,K=1,KMAX) 

WRITE( 92 ) ( ( ( T( I , J ,K) , 1=1 , IMAX ) , J=1 , JMAX) ,K=1 ,KMAX ) 

WRITE(92 )(( (Z(I,J ,K) ,1=1 , IMAX) ,J=1 ,JMAX) ,K=1 , KMAX) 
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K0UHT=K0UNT-1 

WRITE(92)K0UNT,NUM,M, (HSEQ( L) ,L=1 ,4 ) 

write(92)(delx(l) ,l=i,4) ,a,crit,e 

WRITE(92)REV,DLIQ,PSAT,TT,VISC 

WRITE ( 92 ) ( ( (U( I , J ,K) ,1=1 , IMAX) , J=1 , JMAX) ,K=1 ,KMAX) 

WRITE(92)( ( (V(I,J,K) ,1=1, IMAX),J=1, JMAX) ,K=1,KMAX) 

WRITE ( 92 ) ( ( (W( I ,J ,K) ,1=1 ,IMAX) , J=1 , JMAX) ,K=1 ,KMAX) 
WRITE(92)(((P(l,J,K),I=l,IMAX),J=l,JMAX),K=l,KMAX) 

END FILE 92 
REWIND 92 
DO 2 1=1 ,IMAX 
DO 2 J=1 , JMAX 
DO 2 K=1,KMAX 

2 WRE(I,J,K)=SQRT(U(I,J,K)*U(I,J,K)+V(I,J,K)«V(I,J,K)4W(I,J,K)*W(I,J 

i»k) ) 

PRINT 1000 
DO 97 K=1,KMAX 
PRINT 1005 

97 PRINT 1003,((I,J,K,U(I,J,K),V(I,J,K),W(I,J,K),P(I,J,K),WRE(I,J,K), 
IRES ( I , J ,K) , 1=1 , IMAX ) , J=1 , JMAX ) 

GO TO 99 
369 PRINT 1014 
GO TO 98 
371 PRINT 1015 
GO TO 98 
C 

C FORMAT STATEMENTS 
C 

55 F0RMAT(32X,’*## OUTLET AIR ANGLE FOR I=',I2,' CHANGED T0',F10.3) 
700 F0RMAT(9X,I1) 

720 FORMAT (F10 .1 ,F10 . 0,F10 . 5 ,F10 . 5 ,F10 . 0,F10 , 5 ,F10 . 5 ) 

1000 FORMAT (75HO 3 DIMENSIONAL ANALYSIS OF SECONDARY FLOW IN PENN ST 

1ATE 3 BLADE INDUCER ) 

1002 FORMAT ( 8 FI 0 . 7 ) 

1003 FORMAT(3I3,6Fl8.7) 

1004 FORMAT(7I5) 

1005 FORMAT (//2X,1HI,2X,1HJ,2X,1HK,9X,1HU,17X,1HV,17X,1HW,17X,1HP,17X, 
13HVEL,13X,3HRES/) 

1006 FORMAT ( 2X23HAD JUSTMENT SEQUENCE IS ,4ll,38H WHERE U IS 1, V IS 2, 
1 W IS 3, P IS 4.///) 

1007 FORMAT (I3,1HX,I2,1HX,I2,5K GRID,6x,2HM= ,12 ,6x,2HA= ,1PE10 .4//5H RE 
1V= ,1PE10 . 4 , 5X, 5HDLIQ= ,1PE10 . 4 , 5X , 5HPSAT= ,1PE10 . 4 , 5X,3HTT= ,1PE10 . 4 , 
1 5X , 5HVISC= ,1PE10 . 4 / ) 

1012 FORMAT (/6H RELAX, 3X.5HNO OF, 5X.9HTOTAL RMS,7X,7HMAX RMS,7X,7HMIN R 
1MS,16X,44HMAGNITUDE OF BIGGEST ACCEPTED ADJUSTMENT F0R/6H CYCLE, 3X 
1 , 5HTRIES , 3 ( 6X , 8HRESIDUAL ) ,17X ,1HU , 13X ,1HV ,13X ,1HW ,13X ,1HP/ ) 

1013 FORMAT (I6,I8,1P3E14.4,8X,1P4E14.4) 

1014 FORMAT (/20X,9HC0NVERGED/ ) 

1015 FORMAT ( /20X.33HMAXIMUM NUMBER OF CYCLES EXECUTED/) 

C 

C END OF PROGRAM 
C 

99 STOP 
END 
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SUBROUTINE JACOB 
IMPLICIT LOGICAL*! ($) 


CALCULATION OF JACOBIAN COEFFICIENTS FOR GENERAL- TO-CYLINDRICAL 
COORDINATE TRANSFORMATION 

COMMON U(7,7,39),V(7,7,39),W(7,7,39),P(7,7,39),D(7,7,39),R(7,7,39) 
1 ,T(7 ,7,39) ,z(7 ,7,39 ) ,FR,FT,FZ ,DMAX(U ) ,DELX(U ) ,NSEQ( 4 ) ,RES( 8,8,1»0 ) , 
2AR(7,7,39),AT(7,7,39),AZ(7,7,39),BR(7,7,39),BT(7,7,39),BZ(7,7,39), 
3CR(7,7,39) ,CT(7,7,39) ,CZ(7,7,39) .XH.RSTAR, 
4tt,psat,dliq,rev,a,m,nvar,ntr,irstar,dx,i,j,k,ii,jj,kk,imax,jmax,k 
5MAX,VISC,NBD,KLE,KTE,WRE(7,7,39) 

COMMON / $$STRT/ KSTART 
DO 18 K= KSTART ,KMAX 
DO 18 J=l, JMAX 
DO 18 1=1 ,IMAX 
IF (1-1) 3.3.U 

3 RA=R(I+1,J,K)-R(I,J,K) 

TA=T(l+l,J,K)-T(l,J,K) 

ZA=Z(l+l,J,K)-Z(l,J,K) 

GO TO 7 

U IF (i-IMAX) 5,6,6 

5 RA=(R(l+l,J,K)-R(l-l,J,K))/2. 

TA=(T(l+l,J,K)-T(I-l,J,K))/2. 

ZA=(Z(l+l,J,K)-Z(l-l,J,K) )/2. 

GO TO 7 

6 RA=R(I,J,K)-R(I-1,J,K) 

TA=T(I,J,K)-T(I-1,J,K) 

ZA=Z( I,J,K)-Z(l-1 , j,k) 

7 IF (J-1) 8,8,9 

8 RB=R(I,J+1,K)-R(I,J,K) 

TB=T( I , J+l ,K)-T( I , J ,K) 

ZB=Z(I,J+1,K)-Z(I,J,K) 

GO TO 12 

9 IF ( J-JMAX) 10,11,11 

10 RB=(R(l,J+l,K)-R(l,J-l,K))/2. 

TB=(T(l,J+l,K)-T(l,J-l,K))/2. 

ZB=(Z(l,J+l,K)-Z(l,J-l,K))/2. 

GO TO 12 

11 RB=R(I,J,K)-R(I,J-1,K) 

TB=T(l ,J,K)-T( I ,J-1 ,k) 

ZB=Z(I,J,K)-Z(I,J-1,K) 

12 IF (K-l) 13.13.1 1 * 

13 RC=R( I , J ,K+1 )-R( I , J ,K) 

TC=T( I , J ,K+1 )-T( I , J ,K) 

ZC=Z( I ,J,K+1 )-Z( I, J,K) 

GO TO 17 

14 IF (K-KMAX) 15,16,16 

15 RC=(R(l,J,K+l)-R(l,J,K-l))/2. 

TC=(T(l,J,K+l)-T(l,J,K-l))/2. 

ZC=(Z(l,J,K+l)-Z(l,J,K-l))/2. 

GO TO 17 

16 RC=R(I,J,K)-R(I,J,K-1) 

TC=T( I , J ,K)-T( I , J ,K-1 ) 
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ZC=Z( I , J ,K) -Z( I , J ,K-1 ) 

17 B=RA*(TB*ZC-TC*ZB)+RB* (TC*ZA-TA*ZC)+RC* (TA*ZB-TB*ZA) 
AR(l,J,K)=(TB*ZC-TC*ZB)/B 

BR ( I , J ,K) =( TC*ZA-TA*ZC ) /B 
CR(l,J,K)=(TA«ZB-TB*ZA)/B 
AT ( I , J , K ) = ( ZB*RC-ZC*RB ) /B 
BT ( I , J , K) = ( ZC*RA-ZA*RC ) /B 
CT ( I , J ,K) = ( ZA*RB-ZB*RA)/B 
AZU^.KMRB^C-RC^TBj/B 
BZ(I,J,K)=(RC*TA-RA*TC)/B 

18 CZ(l,J,K)=(RA*TB-RB*TA)/B 
RETURN 

END 
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SUBROUTINE RES ID 
IMPLICIT LOGIC AL*1 ($) 

C 

C CALCULATION OF POINT RESIDUAL 
C 

COMMON U(7,7,39),V(7,7,39),W(7,7,39),P(7,7,39),D(7,7,39),R(7,7,39) 
l,T(7,7,39),Z(7,7,39),FR,JT,FZ,DMAX(lt),DELX(l*),NSEQ(l+),RES(8, 8,1*0), 
2AR(7,7,39),AT(7,7,39),AZ(7,7,39),BR(7,7,39),BT(7,7,39),BZ(7,7,39), 
3CR(7 ,7,39) ,CT(7,7,39) ,CZ(7,7,39) .XH.RSTAR, 
i+tt,psat,dliq,rev,a,m,nvar,ntr,irstar,dx,i,j,k,ii,jj,kk,imax,jmax,k 
5MAX,VISC,NBD,KLE,KTE,WRE(7,7,39) 

COMMON /$$LOSS/ TAUR(7,7,26) ,TAUT(7,7,26) ,TAUZ(7,7,26) ,TAU(7,7,26) 
COMMON/ $$CFRE/ SLOPEP( 25 ) ,SLOPES( 25 ) ,BP( 25 ) ,BS ( 25 ) 

IF (NBD) 1+00,1*01,99 
C 

C CHECK WALL BOUNDARY CONDITIONS 
C 

99 IF (K-2) 1*01,100,103 

C 

C K=2 

C 

100 IF ( (l-l)*(lMAX-l) ) U00,121, 12U 

121 IF ( ( J-l )*( JMAX-J ) ) 1*00,102,122 

122 IF (AZ(I,J,K) ) 123,351,123 

123 W(I,J,K)=-U(I,J,K)*AR(I,J,K)/AZ(I,J,K) 

GO TO 351 

12l* IF ( (J-1)«(JMAX-J)) 1*00,111,351 

102 DDD=(BT(I,J,K)/R(I,J,K))/(AR(I,J,K)*BZ(I,J,K)-BR(I,J,K)»AZ(I,J,K)) 
W(I,J,K)=-V(I,J,K)»AR(I,J,K)*DDD 

GO TO 351 

103 IF (K-KMAX) 108,101* ,1*00 
C 

C **« K=KMAX *** 

C 

101* IF ( ( 1-1 ) * ( IMAX-I ) ) 1+00,105,106 

105 u(i,j,k)=-w(i,j,k)*az(i,j,k)/ar(i,j,k) 

106 IF ( ( J-l ) * ( JMAX-J ) ) 1*00,107,351 

107 v(i,j,k)=-r(i,j,k)*(u(i,j,k)»br(i,j,k)+w(i,j,k)*bz(i,j,k))/bt(i,j, 

IK) 

GO TO 351 

108 if(k-kte) 12,8,8 

C 

C '•»» K.GE.KTE 

C 

C FORCE VELOCITIES AND PRESSURES TO SATISFY KUTTA CONDITION 
C 

8 IF((J-1)*(JMAX-J))12,11,12 
11 P(l,J,K)=(P(l,l,K)+P(l,7,K))/2. 

U(l,J,K)=(U(l,l-,K)+U(l,7,K))/2. 

V(l,J,K)=(V(l,l,K)+V(l,7,K))/2. 

W(l,J,K)=(W(l,l,K)+W(l,7,K))/2. 

C 

C K IS NEITHER 2 NOR KMAX *»• 

C 
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12 IP ( (l-l)*(lMAX-l) ) 400,109,110 

109 IF ((J-lWJMAX-J)) 400,112,114 

110 IF ( (J-1)*(JMAX-J) ) 400,111,351 

111 w(i, j , k)=-(u(i ,j ,k)*br(i,j,k)+v(i,j ,k)*bt(i,j ,k)/r(i ,j,k) )/bz(i ,j, 

IK) 

GO TO 351 

112 DDD=(BT(l,J,K)/R(I,J,K))/(AR(l,J,K)*BZ(l,J,K)-BR(l,J,K)*AZ(l,J,K)) 
w(i,j,k)=-v(i,j,k:)«ar(i,j,k)*ddd 
U(l,J,K)=V(l,J,K)«AZ(l,J,K)*DDD 

GO TO 351 

114 U(I, j,k)=-w(i,j,k)*az(i,j,k)/ar(i, J,K) 

351 IF((K.LT.KLE).OR.(K.GT.KTE))GO TO 401 
IF(VISC)402,401 ,402 
C 

C FOR VISCOUS SOLUTION, BLADE SURFACES HAVE ZERO VELOCITY 
C 

402 IF((J-l)*(JMAX-J))400,403,401 

403 U(I,J,K)=0.0 
v(i,j,k)=o.o 
w(i,j,k)=o.o 

c 

C CALCULATE ALL DERIVATIVES 
C 

401 IF( 1-1)400,354, 3 53 

354 UA=(U(l+l , J,K)-U(l ,J,K) ) 

VA=(V(l+l,J ,K)-V(l,J ,K) ) 

WA=(W(I+1,J,K)-W(I,J,K)) 

PA=(P(l+l,J,K)-P(I,J,K)) 

IF (TT) 511,357,511 

511 da=(d(i+i,j,k)-d(i,j ,K) ) 

GO TO 357 

353 IF (1-IMAX) 355,356,400 

356 UA=(U(I,J,K)-U(I-1,J,K)) 

VA= ( V(l , J ,K)-V( 1-1 , J,K) ) 

WA=(W(l,J,K)-W(l-l,J,K) ) 

PA=(P(l,J ,K)-P(l-l,J ,K) ) 

IF (TT) 521,357,521 

521 DA=(D(I,J,K)-D(I-1,J,K)) 

GO TO 357 

355 UA=(U(l+l,J,K)-U(l-l,J,K))/2. 

VA=(V(l+l,J,K)-V(l-l,J,K) )/2. 

WA= ( W( 1+1 , J ,K) -W( 1-1 , J ,K) ) / 2 . 

PA=(P(l+l,J,K)-P(l-l,J,K))/2. 

IF (TT) 531,357,531 

531 DA=(D(l+l,J,K)-D(I-l,J,K))/2. 

357 IF (J-1) 400,359,358 

359 UB=(U(I,J+1,K)-U(I,J,K)) 

VB=(V( I ,J+1 ,K)-V(I ,J,K) ) 

WB=(W(I,J+1,K)-W(I,J,K)) 

pb=(p(i,j+i,k)-p(i,j ,k) ) 

IF (TT) 541,362,541 
541 DB=(D(I,J+1,K)-D(I,J,K)) 

GO TO 362 

358 IF (J-JMAX) 360,361,400 
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361 UB=(U(l,J,K)-U(I,J-l,K)) 

VB=(V(l,J,K)-V(I,J-l,K)) 

WB=(W(l,J,K)-W(l,J-l,K)) 

PB=(P(I ,J,K)-P(I , J-l ,K) ) 

IF (TT) 551,362,551 

551 DB=(D(l,J ,K)-D(I,J-1,K) ) 

GO TO 362 

360 UB=(U(l,J+l,K)-U(l,J-l,K))/2. 

VB=(V(l,J+l,K)-V(l,J-l,K))/2. 

WB=(W(l,J+l,K)-W(l,J-l,K))/2. 

PB= ( P( I , J+l , K) -P( I , J-l , K) ) /2 . 

IF (TT) 561,362,561 
561 DB=(D(l,J+l,K)-D(l,J-l,K))/2. 

362 IF (K-l) 1*00, 36U, 363 

361* UC=(U(l ,J,K+1 )-U(l , J,K) ) 

VC=(V(I,J,K+1)-V(I,J,K)) 

WC=(W(l,J,K+l)-W(I,J,K)) 

PC=(P( I , J,K+1 )-F( I ,J,K) ) 

IF (TT) 571,367,571 
571 dc=(d(i,j,k+i)-d(i,j,k) ) 

GO TO 367 

363 IF (K-KMAX) 365,366,1*00 

366 UC=(U(l ,J,K)-U( I ,J,K-1 ) ) 

VC=(V(I,J,K)-V(I,J,K-1)) 

WC=(W(I,J,K)-W(I,J,K-1)) 

PC=(P(I,J,K)-P(I,J,K-1)) 

IF (TT) 581,367,581 

581 DC=(D(l ,J,K)-D(l ,J,K-1 ) ) 

GO TO 367 

365 UC=(U(l,J,K+l)-U(l,J,K-l))/2. 

VC=(V(l,J,K+l)-V(l,J,K-l))/2. 

WC= ( W( I , J ,K+1 ) -W( I , J ,K-1 ) ) /2 . 

PC=(P(l,J,K+l)-P(l,J,K-l))/2. 

IF (TT) 591,367,591 
591 DC=(D(l,J,K+l)-D(l,J,K-l))/2. 

C 

C CONVERT ALL DERIVATIVES FROM GENERAL TO CYLINDRICAL COORDINATES 
C 

367 UR=AR(I,J,K)*UA+BR(I,J,K)*UB+CR(I,J,K)*UC 
UT=AT(I,J,K)*UA+BT(I,J,K)*UB+CT(I,J,K)»UC 
UZ=AZ(I,J,K)*UA+BZ(I,J,K)*UB+CZ(I,J,K)«UC 
VR=AR( I , J ,K) *VA+BR (l,J,K) #VB+CR( I ,J,K)*VC 
VT=AT(l,J,K)*VA+BT(l,J,K)*VB+CT(l,J,K)*VC 
VZ=AZ ( I , J , K ) *VA+BZ ( I , J , K) *VB+CZ ( I , J ,K) *VC 
WR=AR( I , J ,K) *WA+BR ( I , J ,K)*WB+CR(l , J ,K)*WC 
WT=AT(l,J,K)*WA+BT(l,J,K)»WB+CT(l,J,K)*WC 
WZ=AZ( I , J ,K) *WA+BZ(l , J ,K)»WB+CZ( I , J ,K)«WC 
PR=AR(I,J,K)*PA+BR(I,J,K)*PB+CR(I,J,K)*PC 
PT=AT(l,J,K)*PA+BT(l,J,K)*PB+CT(l,J,K)*PC 
PZ=AZ(l,J,K)»PA+BZ(l,J,K)*PB+CZ(l,J,K)«PC 
IF (TT) 370,375,370 

370 DR=AR(I,J,K)*DA+BR(I,J,K)*DB+CR(I,J,K)*DC 
DT=AT(l ,J ,K)»DA+BT(I ,J ,K)*DB+CT(l ,J,K)*DC 
DZ=AZ(l,J,K)*DA+BZ(l,J,K)*DB+CZ(l,J,K)«DC 
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GO TO 380 
375 DR=0.0 
DT=0.0 
DZ=0.0 

380 UU=U(I,J,K) 

W=V(l,J,K) 

WW=W(l,J,K) 

DD=D(l,J,K) 

RR=R(l,J,K) 

CALCULATE LOSS TERMS 

IF (VISC )385, 390, 385 

385 if( (k.lt.kle) .or. (k.gt.kte) )go to 390 

CALL DLOSS 
IF (I-1)U00,501,502 

501 TRA=TAUR(I+1,J,K)-TAUR(I,J,K) 

TTA=TAUT ( I +1 , J , K ) -TAUT ( I , J , K ) 

TZA=TAUZ ( 1+1 , J ,K)-TAUZ(I ,J ,K) 

GO TO 503 

502 IF(l-IMAX)50l+, 505, 1+00 

505 TRA=TAUR(I,J,K)-TAUR(I-1,J,K) 
TTA=TAUT(I , J,K)-TAUT( 1-1 , J,K) 
TZA=TAUZ( I , J ,K) -TAUZ( 1-1 , J ,K) 

GO TO 503 

50l+ TRA=(TAUR(I+1,J,K)-TAUR(I-1,J,K) )/2. 
TTA= ( TAUT ( 1+1 , J ,K)-TAUT( 1-1 ,J ,K) )/2. 
TZA=(TAUZ(l+l,J,K)-TAUZ(I-l,J,K))/2. 

503 IF(j-l)l+00,506,507 

506 TRB=TAUR(I,J+1,K)-TAUR(I,J,K) 

itb=taut(i,j+i,k)-taut(i,j ,k) 

TZB=TAUZ(I,J+1,K)-TAUZ(I,J,K) 

GO TO 508 

507 IF(J-JMAX)509, 510, 1+00 

510 TRB=TAUR ( I , J , K ) -TAUR ( I , J-l , K ) 
TTB=TAUT( I , J ,K) -TAUT( I , J-l ,K) 

TZB=TAUZ ( I , J , K) -TAUZ ( I , J-l ,K) 

GO TO 508 

509 TRB=(TAUR(I,J+l,K)-TAUR(I,J-l,K))/2. 
TTB= ( TAUT ( I , J+l , K ) -TAUT ( I , J-l , K) ) / 2 . 
TZB= ( TAUZ ( I , J+l ,K)-TAUZ( I , J-l ,K) )/2 . 

508 IF(K-KLE)390, 518 ,512 

518 TRC=TAUR(I,J,K+1)-TAUR(I,J,K) 
TTC=TAUT ( I , J , K+l ) -TAUT ( I , J , K ) 
TZC=TAUZ ( I , J ,K+1 ) -TAUZ ( I , J ,K) 

GO TO 513 

512 IF(K-KTE) 51^ >515 ,390 
515 trc=taur(i,j,k)-taur(i,j,k-i) 
TTC=TAUT(I , J ,K) -TAUT ( I , J ,K-1 ) 
TZC=TAUZ( I , J ,K)-TAUZ( I , J ,K-1 ) 

GO TO 513 

511+ TRC= ( TAUR ( I , J ,K+1 ) -TAUR (I ,J ,K-1 ) )/2 . 
TTC=(TAUT(I , J ,K+1 )-TAUT( I , J ,K-1 ) )/2 . 
TZC=(TAUZ( I , J ,K+1 )-TAUZ( I ,J ,K-1 ) )/2 . 
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513 TRWZ=AZ(l,J,K)*TRA+BZ(l,J,K)«TRB+CZ(l,J,K)»TRC 
TTWZ=AZ( I , J ,K) *TTA+BZ( I , J ,K) «TTB+CZ ( I , J ,K)*TTC 
TZWTsATUjJ.K^TZA+BTte.J.Kj^TZB+CTCl.J.K^TZC 
FR—TRWZ/DD 
FT=-TTWZ/DD 
FZ=-TZWT/(DD*RR) 

GO TO 391 

390 FR=0. 0 
FT=0.0 
FZ=0 . 0 

CALCULATE POINT RESIDUALS 

391 Rl=PR/DD+UU*UR+VV*UT/RR+WW*UZ-( ( W+RR*REV )* ( W+RR*REV ) ) /RR+FR 
R2=PT/ ( DD*RR ) +UU*VR+W»VT/RR+WW*VZ+UU*W/RR+2 . *UU*REV+FT 
R3=PZ/DD+UU*WR+W*WT/RR+WW*WZ+FZ 

Rl+ =UU/ RR+UR+VT/ RR+WZ+ ( UU*DR+W*DT/RR+WW«DZ ) /DD 
RES ( I , J , K) =R1*R1+R2*R2+R3*R3+R4 *R4 
400 RETURN 
END 
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SUBROUTINE DLOSS 
IMPLICIT LOGICAL*! ($) 
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CALCULATION OF LOSS TERMS IN MOMENTUM EQUATIONS 

COMMON U(7,7,39),V(7,7,39),W(7,7,39),P(7,7,39),D(7,7,39),R(7,T,39) 
1 ,T( 7 ,7 , 39 ) ,Z( 7 ,7 , 39 ) ,FR,FT,FZ ,DMAX( U) ,DELX( 10 ,NSEQ( 1+ ) ,RES ( 8 ,8 ,1*0 ) , 
2AR(7,7,39),AT(7,7,39),AZ(7,7,39),BR(7,7,39),BT(7,7,39),BZ(7,7,39), 
3CR(7,7 ,39) ,CT(7 ,7 ,39) >CZ(7 ,7 »39) ,XH,RSTAR, 

1+TT,PSAT,DLIQ,REV,A,M,NVAR,NTR,IRSTAR,DX,I,J,K,II,JJ,KK,IMAX,JMAX,K 

5MAX,VISC,NBD,KLE,KTE,WRE(7,7,39) 

common /$$loss/ taur( 7 , 7 , 26 ) ,taut( 7 , 7 , 26 ) ,tauz( 7 , 7 , 26 ) ,tau( 7 , 7 , 26 ) 

COMMON/ $$CFRE/ SL0PEP( 25) .SLOPES (25) ,BP(25) ,BS(25) 

DIMENSION VEL(7) 

OMEGA AND RTIP ARE SPECIFIC PARAMETERS OF PENN STATE INDUCER 

0MEGA=1+ 50 . *3 .11*159/30 . 

RTIP=18. 25/12. 

DEN=D(I,J,K) 

JMIN=1 
UBAR=0 . 0 
UTIP=RTIP*OMEGA 

FIND PASSAGE AVERAGED VELOCITY 

DO 100 Jl=l , JMAX 

VEL(J1)=SQRT(U(I,J1,K)*U(I,J1,K)+V(I,J1,K)*V(I,J1,K)+W(I,J1,K)*W(I 
1 ,J1 ,K) ) 

UBAR=UBAR+VEL( J1 ) 

100 CONTINUE 
XJMAX=JMAX 
UBAR=UBAR/ XJMAX 

CALCULATE REYNOLDS NUMBER, FIND APPROPRIATE SKIN FRICTION 

COEFFICIENT, THEN CALCULATE THE COMPONENTS OF SHEAR STRESS 

UU=UBAR*UTIP 

rr=r(i,j,k)*rtip 

RE=UU*RR/VISC 

CF1=BP ( K ) / ( RE**SLOPEP ( K ) ) 

CF2=BS ( K ) / ( RE**SLOPES ( K ) ) 

TAU1=CF1*DEN*UBAR**2 . / 2 . 

TAU2=-CF2*DEN*UBAR**2 . / 2 . 

XNUM=J-JMIN 
XDEN=JMAX-JMIN 
RATIO=XNUM/ XDEN 

TAU( I , J ,K)=TAU1-RATI0*(TAU1-TAU2 ) 

IF(J-l) 1+00, 359, 358 
359 DELU=(U( I , J+l ,K)-U( I , J ,K) ) 

DELV=(V( I ,J+1 ,K)-V( I , J,K) ) 

DELW=(W(l ,J+1 ,K)-W(l ,J,K) ) 

DELVEL= ( VEL (J+l )-VEL( J ) ) 

GO TO 362 
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358 IF ( J-JMAX) 360,361,400 

361 DELU=(U(l,J,K)-U(l,J-l,K)) 
DELV=(V(l,J,K)-V(l,J-l,K)) 

DELW=(W( I , J ,K)-W( I , J-l ,K) ) 

DELVEL= ( VEL ( J ) - VEL ( J-l ) ) 

GO TO 362 

360 DELU= (U( I , J+l ,K) -U( I , J-l ,K) )/2. 
DELV=(V(l,J+l,K)-V(l,J-l,K))/2. 
DELW=(W(l,J+l,K)-W(l,J-l,K))/2. 
DELVEL= ( VFL( J+l )-VEL( J-l ) ) /2 . 

362 TAUR( I , J ,K)=TAU( I , J , K ) *DELU/ DELVEL 
TAUT ( I ,J,K)=TAU(I,J,K)*DELV/DELVEL 
TAUZ ( I , J , K) =TAU( I , J , K ) *DELW/ DEL VEL 

bOO RETURN 
END 



SUBROUTINE ADJ 
IMPLICIT LOGICAL*! ($) 
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C 

C REDUCTION OF TOTAL RESIDUAL BY SUCCESSIVE VARIATIONS OF 
C PARAMETERS U, V, W, P 

C 

COMMON U(7,7,39),V(7,7,39),W(7,7,39),P(7,7,39),D(7,7,39),R(7,7,39) 
1,T(7 ,7 *39) ,Z(7 ,7 ,39) ,FR,FT,FZ,DMAX(4) ,DELX(4) ,NSEQ(4) ,RES(8,8,40) , 
2AR(7,7,39),AT(7,7,39),AZ(7,7,39),BR(7,7,39),BT(7,7,39),BZ(7,7,39), 
3CR(7,7,39) ,CT(7,7,39) ,CZ(7,7,39) .XH.RSTAR, 

4TT,PSAT,DLIQ,REV,A,M,NVAR,NTR,IRSTAR,DX,I,J,K,II,JJ,KK,IMAX,JMAX,K 

5MAX,VISC,NBD,KLE,KTE,WRE(7,7,39) 

C 

C CURRENT VALUES OF ALL AFFECTED QUANTITIES TEMPORARILY STORED 
C 

Hl=RES(lI,JJ,KK) 

IF (II-1) 1101,1101,1100 

1100 H2=RES(II-1,JJ,KK) 

1101 IF (II-IMAX) 1102,1103,1103 

1102 II3=RES(II+1,JJ,KK) 

1103 IF (JJ-1) 1105, 1105, 110U 

1104 H4=RES(II,JJ-1,KK) 

1105 IF ( JJ-JMAX) 1106,1107,1107 

1106 H5=RFS(II,JJ+1,KK) 

1107 IF (KK-1 ) 1109,1109,1108 

1108 H6=RES(II,JJ,KK-1) 

1109 IF (KK-ro-tAX) 1110,1111,1111 

1110 H7=RES(II,JJ,KK+1) 

1111 RHLD=RSTAR 
DX=DELX ( NVAR ) 

HU=U(II,JJ,KK) 

HV=V( II ,JJ,KK) 

HW=W(lI ,JJ ,KK) 

HP=P( II , JJ ,KK) 

HD=D(II,JJ,KK) 

C 

C SUCCESSIVELY APPLY TRIAL VARIATIONS TO U, V, W, P 
C 

5 DO 480 MA=1,M 

420 GO TO (422, 423, 424, 425), NVAR 

422 U(II,JJ,KK)=HU+DX 
GO TO 421 

423 V(II,JJ,KK)=HV+DX 
GO TO 421 

424 W( II ,JJ,KK)=HW+DX 
GO TO 421 

425 P( II ,JJ,KK)=HP+DX 
IF (TT) 461,421,461 

4 6l I=II 
J=JJ 
K=KK 

CALL STATE 

421 NTR=NTR+1 
IRSTAR=2 


o o o 


158 


CALL STAR 

IF (RHLD-RSTAR) 1*30,430,431 

431 IF (ABS(DX)-DMAX(NVAR)) 457,457,920 
920 DMAX(NVAR)=ABS(DX) 

GO TO 457 

430 IF (DX) 433,457,432 

432 DX=-DX 

GO TO 420 

433 DX=— A*DX 
480 CONTINUE 

RESTORE ALL AFFECTED QUANTITIES TO ORIGINAL VALUES 

U(II,JJ,KK)=HU 

V(II,JJ,KK)=HV 

W(lI,JJ,KK)=HW 

P(II,JJ,KK)=HP 

D( II , JJ ,kk)=hd 

445 RSTAR=RHLD 

RES ( II , JJ,KK)=H1 

IF (Il-l) 446,446,447 

447 RES(II-1,JJ,KK)=H2 

446 IF (II-IMAX) 449,448,448 

449 RES(II+1,JJ,KK)=H3 

448 IF (JJ-1) 450,450,451 

451 RES(II,JJ-1,KK)=h4 

450 IF (JJ-JMAX) 453,452,452 

453 RES ( II ,JJ+1 ,KK)=H5 

452 IF (KFC-1 ) 454,454,455 

455 RES(II,JJ,KK-1)=H6 

454 IF (KK-KMAX) 456,457,457 

456 RES ( II ,JJ,KK+1 )=H7 

457 RETURN 
END 
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SUBROUTINE STAR 
IMPLICIT L0GICAL*1 ($) 

CALCULATION OF LOCAL STAR RESIDUAL 

COMMON U(7 ,7 ,39) ,V(7 ,7 ,39) ,W(7,7,39) ,P(7,7,39) ,D(7,7,39) ,R(7,7,39) 
1»T(7 ,7*39) ,z( 7 »7 ,39) ,FR,FT,FZ,DMAX(.4) ,DELX(4) ,NSEQ(4) ,RES(8,8,40) , 
2AR(7,7,39) ,AT(7,7,39) ,AZ(7,7,39) ,BR(7,7,39) ,BT(7,7,39) ,BZ(7,7,39) , 
3CR(7,7»39) ,CT(7»7»39) ,CZ(7,7,39) ,XH,RSTAR, 

4TT,PSAT t DLIQ,REV,A > M,NVAR,NTR,IRSTAR,DX > I,J,K,II,JJ,KK,IMAX,JMAX,K 
5MAX,VISC ,NBD ,KLE,KTE,WRE(7*7 ,39) 

RSTAR=0 . 0 

IG0=0 

I=II 

J=JJ 

K=KK 

NBD=1 

389 GO TO (375,378) ,IRSTAR 

378 CALL RES ID 

375 RSTAR=RSTAR+RES( I , J ,K) 

NBD=0 

400 IGO=IGO+l 

GO TO (391,392,393,394,395,396, 402),IG0 

391 1=11-1 

IF (1-1) 400,389,389 

392 1=11+1 

IF (i-IMAX) 389,389,400 

393 1=11 
J=JJ-1 

IF (J-1) 400,389,389 

394 J=JJ+1 

IF (J-JMAX) 389,389,400 

395 J=JJ 
K=KK-1 

IF (K-l) 400,389,389 

396 K=KK+1 

IF (K-KMAX) 389,389,400 
402 RETURN 
END 



160 


SUBROUTINE STATE 
IMPLICIT LOGICAL*! ($) 


C 

C CALCULATION OF POINT DENSITY 
C 

COMMON U(7,7,39) ,v(7,7,39) ,W(7,7,39) ,P(7,7,39) ,D(7,7,39) ,R(7,7,39) 
l,T(7,7,39),Z(7,7,39),FR,FT,FZ,DMAX(lt),DELX(U),NSEQ(U),RES(8,8,U0), 
2AR(7,7,39) »AT(7,7 *39) »AZ(7 *7 *39) ,BR(7,7,39) ,BT(7,7,39) ,BZ(7,7,39) , 
3CR(7,7,39) »CT(7 ,7 ,39) »CZ(7 »7 »39) ,XH,RSTAR, 

4TT,PSAT,DLIQ,REV,A,M,NVAR,NTR,IRSTAR,DX,I,J,K,II,JJ,KK,IMAX,JMAX,K 
5MAX , VISC ,NBD ,KLE ,KTE ,WRE( 7 ,7 , 39 ) 

PSP=PSAT-P( I , J ,K) 

IF (PSP) 101,101,102 

101 D(I,J,K)=DLIQ 
GO TO 103 

102 D(I,J,K)=DLIQ/(1.+TT*PSP) 

103 RETURN 
END 
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APPENDIX C 


Fortran Listing of Program to Solve Hotwire Equations 


ROTATING HOTWIRE MEASUREMENTS IN PENN STATE INDUCER PALL 1973 

PROGRAM DEVELOPED TO CALCULATE RELATIVE MEAN VELOCITIES, 
TURBULENCE INTENSITIES AND STRESSES FROM TOE EXPERIMENTAL DATA 

IMPLICIT L0GICAL*1 ($) 

DIMENSION rV0LTl(U,6,10),FVOLT2(U,6,10),FVOLT3(M,10) 

DIMENSION SL0PE1( 1* ) ,SL0PE2( U ) ,SL0PE3(U ) 

DIMENSION FVELl(6,10),FVEL2(6,10),FVEL3(6,10) 

DIMENSION Vl(lt,15),V2(U,15),V3(U,15) 

DIMENSION NPTS(15),CINCH(1*,15),CALV1(1(,15),CALV2(1*,15),CVEL(M5) 

DIMENSION VCHAN2(l*,6,10),VELl(6,10),VEL2(6,10) 

DIMENSION VCHAN3(*t,6,10),VEL3(6,10),CALE3(U,15),CALV3(l*,15) 
DIMENSION CALE1(U,15),CALE2(1*,15),RRAD(6) 

DIMENSION THETA(IO) ,RAD(6) ,VCHAN1(1*,6,10) 

DIMENSION C(3»3) ,A(6,6) ,TEMP(10) 

REAL*8 CC(3,3),AA(6,6),D(3),B(6),DET 

DEFINE PROGRAM CONSTANTS 

UTIP=1*50 . *2 . *3 . 1»*159*18 . 25/ ( 60 . *12 . ) 

FACT=SQRT(2.)/2. 

FACT0R=19./B. 

ATTEN=27.1* 

RAT=3. 1^159/180. 

READ IN CONVERGENCE CRITERIA 

READ ( 5 , 82 ) NREPS ,EPS 

READ IN HOTWIRE CONSTAITTS FOR L/D 

READ ( 5 , 70 ) XEL , XK2 , XK3 

READ IN CALIBRATION CURVES 

READ ( 5 , 6 ) CALA1 , CALA2 , CALA3 

CALA1=CALA1*RAT 

CALA2«CALA2*RAT 

CALA3=CALA3*RAT 

READ(5,5)NCAL 

DO 100 I-l.NCAL 

READ( 5 ,10 )NPTS ( I ) ,TEMP( I ) 

JJ=NPTS(I) 

DO 105 J=1,JJ 

READ(5,15)CINCH(I,J),CALE1(I,J),CALE2(I,J),CALE3(I,J) 

CVEL( I , J)*66 . 7*SQRT( CINCH ( I , J) ) 

APPLY TEMPERATURE CORRECTION TO CALIBRATION CURVES 

CVEL(I,J)-CVEL(I,J)*(1.0+0.001*(TEMP(.T)-73.0)) 

CALEl(I,J)«CALEl(I,J)*(l.O+0.00ll**(TEMP(j)-73.0)) 

CALE2 ( I , J ) -CALE2 ( I , J ) * ( 1 . 0+0 . 0011** ( TEMP( J ) -73 . 0 ) ) 
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CALE3 ( I , j ) =CALE3 ( I , J ) * ( 1 . 0+0 . OOlM* ( TEMP ( J ) -73 . 0 ) ) 

CORRECT CALIBRATION CURVE FOR SKEWED PROBE 

VI ( I , J )=SQRT ( ( CVEL( I , J ) *C0S ( CALA1 ) ) »*2 . +XKL*XK1« ( CVEL( I , J ) »SIN ( CAL 
1A1) )**2. ) 

V2 ( I , J ) =SQRT( ( CVEL ( I , J ) *C0S ( CALA2 ) ) ** 2 . +XK2«XK2* ( CVEL( I , J ) «SIN ( CAL 
1A2 ) ) **2 . ) 

V3(I,J)=SQRT( ( CVEL ( I , J ) *C0S ( CALA3 ) ) **2 . +XK3*XK3* ( CVEL ( I , J ) *SIN ( CAL 
1A3))*«2.) 

105 CONTINUE 

100 CONTINUE 

USE VOLTAGE SQUARED VS. ROOT VELOCITY CURVE FOR CALCULATIONS 

CALV1 ( I , J ) =CALE1 ( I , J ) **2 . 

CALV2 ( I , J ) =CALE2 ( I , J ) **2 . 

CALV3( I ,J)=CALE3(I ,J)**2. 

Vl(l,J)=SQRT(Vl(I,J)) 

V2(I,J)=3QRT(V2(l,J)) 

V3(l,J)=SQRT(V3(l,J) ) 

102 CONTINUE 

101 CONTINUE 

DO 103 1=1 ,NCAL 

SL0PEl(l)=(CAI.Vl(I,U)-CALVl(I,2))/(Vl(l,lt)-Vl(I,2)) 

SL0PE2(l) = (CALV2(l,lt)-CALV2(l,2))/(V2(l,lt)-V2(I,?)) 
SL0PE3(I)=(CALV3(I,^)-CALV3(I,2))/(V3(I,U)-V3(I,2)) 

103 CONTINUE 

SPECIFY RADIAL AND TANGENTIAL MEASUREMENT LOCIATIONS 
READ(5,20)N0R,N0T 

READ( 5 , 25 ) ( THETA( I ) , 1=1 , NOT) 

REAU( 5 » 32 ) (RAD(l) ,1=1, NOR) 

DO 112 1=1, NOR 

RRAD(I)=(l8.25-RAD(l))/l8.25 
112 CONTINUE 

READ IN DIRECTION COSINES - THESE CAN BE CALCULATED WITHIN THE 
PROGRAM IF CODED TO DO SO 

READ ( 5 » 70 )A1 ,B1 , Cl 

READ(5,70)A2,B2,C2 

READ(5,70)A3,B3,C3 

READ(5,70)D1,E1,F1 

READ(5,70)D2,E2,F2 

READ(5,70)D3,E3,E3 

READ ( 5 , 7 0 ) G1 , XH1 , XII 

READ ( 5 , 70 ) G2 ,XH2 , XI2 

READ(5,70)G3,XH3,XI3 


CALCULATE EQUATION COEFFICIENTS 

AU=XKL»XKL*A1*A1+A2*A2+A3*A3 
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B4=XK1»XK1»B1»B1+B2»B2+B3*B3 
C4=XK1*XK1*C1»C1+C2*C2+C3*C3 
D4=2 . * ( XKL*XK1*A1*B1+A2*B2+A3*B3 ) 

E4=2 . * ( XK1*XK1*A1*C1+A2»C2+A3»C3 ) 

F4=2 . * ( Xja*xm*Bl*Cl+B2*C2+B3*C3 ) 

A5=XK2«XK2«D1»D1+D2»D2+D3*D3 
B5=XK2*XK2»E1»E1+E2*E2+E3*E3 
C5=XK2»XK2*F1*F1+F2*F2+F3*F3 
D5=2 • ( XK2*XK2*D1»E1+D2*E2+D3*E3 ) 

E5=2.*(XK2«XK2»D1*F1+D2»F2+D3»F3) 

F5=2 . * ( XK2*XK2«E1*F1+E2»F2+E3*F3 ) 

A6=XK3*XK3*G1«G1+G2*G2+G3*G3 
b6=XK3*XK3*XH1*XH1+XH2*XH2+XH3*XH3 
C6=XK3*XK3*XI1»XI1+XI2*XI2+XI3*XI3 
d6=2 . * ( XK3*XK3*G1*XH1+G2*XH2+G3*XII3 ) 

E6=2 . * ( XK3*XK3*G1*XI1+G2*XI2+G3*XI3 ) 

F6=2 . * ( XK3*XK3*XH1*XI1+XH2»XI2+XH3*XI3) 

C(l,l)=SQRT(A4) 

C(l,2)-D4/(2.»C(l,l)) 

C(1,3)=e4/(2.*C(1,1)) 

C(2,l)=SQRT(A5) 

C(2,2)=D5/(2.*C(2,1)) 

C(2,3)=E5/(2.»C(2,1)) 

C(3,l)=SQRT(A6) 

C(3,2)=D6/(2.*C(3,l) ) 

C(3,3)=E6/(2.*C(3,1)) 

READ IN MEAN AND FLUCTUATING VOLTAGES, FIND CORRESPONDING VELOCITIES 

DO 115 J=l,NOT 
READ ( 5 , 4 0 ) IPRO , JRAD 
DO 110 1=1 ,JRAD 
I F ( IPRO . EQ. . 9 ) GO TO 140 
IPROBE=IPRO 

READ (5,46 ) VCH1 , VCH2 ,VCH3 .EXTEMP 

READ (5,45) FV1SQ , FV2SQ , FV3SQ 

READ( 5,75) FV1P2 , FV1M2 ,FV1P3 ,FV1M3 ,FV2P3 ,FV2M3 

GO TO 145 

14 0 READ( 5 , 47) VCII1 ,VCK2 ,VCH3 , EXTEMP, IPR03E 
READ ( 5 , 4 5 ) FV1SQ , FV2SQ , FV3SQ 

READ (5,75) FV1P2 , FV1M2 , FV1P3 , FV1M3 , FV2P3 , FV2M3 
1 •’ ; 5 VCIIAN1 ( IPROBE , I , J ) =VCH1 
VCHAN2 ( IPROBE , I , J ) =VCH2 
VCHAN3 ( IPROBE , I ,J ) =VCH3 

APPLY TEMPERATURE CORRECTION TO EXPERIMENTAL DATA 
TCORR=l. 0+0. 0028* (EXTEMP-7 3.0) 

VCHAN1 ( IPROBE , I , J ) =TC0RR*VCIIA!!1 ( IPROBE , I , J ) **2 . 

VCHAN2( IPROBE , I , J ) =TC0RR*VCHAN2 ( IPROBE , I ,J)**2. 

VCHAN3 ( IPROBE , I , J ) =TC0RR*VCHAN3 ( IPROBE ,I,J)**2. 
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J J=NPTR ( IPROBE ) 

DO 120 K=1,JJ 

IF(VC!IANl(lPROBE,I,J) .GT.CALVl(IPROBE,K) )G0 TO 120 
RATI 0= ( VCHAIT1 ( IPROBE , I , J)-CALV1 ( IPROBE ,K-1 ) ) / ( CALV1 ( IPROBE ,K) -CALV 
11 ( IPROBE ,K-1 ) ) 

VELl(l,J) =V1 ( IPROBE , K-l ) +RATI 0* ( VI ( IPROBE , K ) -VI ( IPR03E , K-l ) ) 
VELl(l,J)=VELl(l,J)**2. 

K1=K 

GO TO 121 

120 CONTINUE 

RATIO= ( VCKAN1 ( IPROBE , I , J ) -CALV1 ( IPROBE, JJ-1 ) ) / ( CALV1 ( IPROBE , J J ) -CA 
1LV1 ( IPROBE, JJ-1)) 

VEL1 ( I , J ) =V1 ( IPROBE , JJ-1 ) +RATIO* ( VI ( IPROBE , JJ ) -VI ( IPROBE , JJ-1 ) ) 
VELl(l,J)=VELl(l,j)**2. 

Kl=99 

121 JJ=NPTS( IPROBE) 

DO 125 K=1,JJ 

IF ( VCEA7T2 (IPROBE, I, J).GT.CALV2( IPROBE ,K)) GO TO 125 
RATIO= ( VCHAN2 ( IPR03E , I , J ) -CALV2 ( IPROBE ,K-1 ) ) / ( CALV2 ( IPROBE ,K) -CALV 
12 (IPROBE, K-l) ) 

VEL2 ( I , J ) =V2 ( IPROBE , K-l ) +RATIO* ( V2 ( IPROBE , K ) -V2 ( IPROBE , K-l ) ) 
VEL2(I,J)=VEL2(I,J)**2. 

K2=K 

GO TO 122 
125 CONTINUE 

RATIO= ( VCIIAN2 ( IPROBE , I , J ) -CALV2 ( IPROBE , JJ-1 ) ) / ( CALV2 ( IPROBE , JJ ) -CA 
1LV2( IPROBE, JJ-1)) 

VEL2 ( I , J ) =V2 ( IPROBE , JJ-1 ) +RATI 0* ( V2 ( IPROBE , JJ ) -V2 ( IPROBE , JJ-1 ) ) 
VEL2(I,J)=VEL2(I,J)**2. 

K2=99 

122 JJ=NPTS( IPROBE) 

DO 130 K=l, JJ 

IF ( VCHA1I3 (IPROBE, I, J).GT.CALV3( IPROBE ,K)) GO TO 130 
RATI0=(VCHA!T3( IPROBE , I , J ) -CALV3 ( IPROBE , K-l ) ) /(CALV3( IPROBE ,K) -CALV 
13(IPR0BE,K-1)) 

VEL3 ( I , J ) =V3 ( IPROBE , K-l ) +RATIO*( V3 ( IPROBE , K ) -V3 ( IPROBE , K-l ) ) 
VEL3(I,J)=VEL3(I,J)**2. 

K3=K 

GO TO 135 
130 CONTINUE 

RATI 0= ( VCIIAN3 ( IPROBE , I , J ) -CALV3 ( IPROBE , JJ-1 ) ) / ( CALV 3 ( IPROBE , JJ ) -CA 
1LV3( IPROBE, JJ-1)) 

VEL3 ( I , J ) =V3 ( IPROBE , JJ-1 ) +RATI 0* ( V3 ( IPROBE , J J ) -V3 ( IPROBE , JJ-1 ) ) 
VEL3(I,J)=VEL3(I,J)**2. 

K3=99 


USE ITEVJTON-RAPHS ON METHOD TO CALCULATE MEAN VELOCITIES 


FIND INITIAL VALtJES OF VELOCITY USING LINEAR EQUATIONS 
C 

135 D(1)=VEL1(I,J) 

D(2)=VEL2(I,J) 

D(3)=VEL3(I,J) 

C 
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C RE-INITIALIZE C(l,J) 

C 

C(l,l)=SQRT(A4) 

C(1,2)=D4/(2.»C(1,D) 

C(1,3)=E4/(2.«C(1,1)) 

C(2,1)=SQRT(A5) 

C(2,2)=D5/(2.*C(2,1)) 

C(2,3)=E5/(2.*C(2,1)) 

C(3,1)=SQRT(A6) 

C(3,2)=D6/(2.*C(3,1)) 

C(3,3)=E6/(2.*C(3,1)) 

DO 111 II®1,3 
DO 111 JJ=1,3 
CC(II,JJ)=C(II,JJ) 

111 CONTINUE 

CALL DLEQD(CC,D,3,1,3,3,DET) 

U0=D(2) 

V0=D(1) 

W0=D(3) 

VINIT=VO 

UINIT=UO 

WINIT=W0 

C 

C BEGIN NEWTON-RAPHSON ITERATION ON NONLINEAR EQUATIONS 
C 

DO 200 N=1,NREPS 

aai=sqrt(a4) 

AA2=SQRT(A5) 

AA3=SQRT(a6) 

BB1=D4 / ( 2 . * AA1 ) 

BB2=D5/(2.«AA2) 

BB3=D6/(2.»AA3) 

CC1=AA1* ( B4/ ( 2 . *A4 )-D4*D4/ ( 8 . *A4*A4 ) ) 

DD1=AA1* ( F4/ ( 2 . *Ak ) -D4*E4 / ( 4 . *A4*A4 ) ) 

EE1=E4 / ( 2 . *AA1 ) 

FF1=AA1* ( C4 / ( 2 . * A4 ) -E4 *E4 / ( 8 . *A4 *A4 ) ) 

CC2=AA2* ( B5 / ( 2 . *A5 ) -D5*D5 / ( 8 . *A5*A5 ) ) 

DD2=AA2* ( F5/ ( 2 . *A5 ) -D5*E5 / ( 4 . *A5*A5 ) ) 

EE2=E5/(2.»AA2) 

FF2=AA2»(C5/(2.»A5)-E5*E5/(8.*A5*A5)) 

CC3=AA3» ( B6/ ( 2 . »A6 ) -D6*D6/ ( 8 . »A6«A6 ) ) 

DD3=AA3* ( F6/(2.*A6) -d6*e6 / ( 4 . *a6*A6 ) ) 

EE3=e6/(2.»AA3) 

FF3=AA3* (c6/(2.*A6)-E6»e6/(8. *A6*A 6 ) ) 

GG1=AA1*V0+BB1*U0+CC1*U0*U0/V0+DD1*U0*W0/V0+FE1*W0+FF1«W0*W0/V0 

l-VELl(l,J) 

GG2=AA2*V0+BB2*U0+CC2*U0*U0/V0+DD2*U0*W0/V0+EE2*W0+FF2*W0*W0/V0 

1-VEL2(I,J) 

GG3=AA3»V0+BB3*U0+CC3*U0*U0/V0+DD3*U0*W0/V0+EE3*W0+FF3*W0»W0/V0 

1-VEL3(I,J) 

G1WV=AA1-CC1*U0*U0/ ( VO*VO ) -DD1*U0*W0/ (VO* VO ) -FF1*W0*W0/ ( VO*VO ) 
G2WV=AA2-CC2*U0*U0/ ( VO*VO ) -DD2*U0*W0/ ( VO*VO ) -FF2*W0*W0/ ( VO*VO ) 
G3WV=AA3-CC3*U0*U0/ ( V0*V0 ) -DD3*U0*W0/ ( V0*V0 ) -FF3*W0*W0/ ( V0*V0 ) 
GlWU=BBl+2 . *CC1*U0/V0+DD1«W0/V0 
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G2WU=BB2+2 . *CC2*UO/VO+DD2«WO/VO 

G3WU=BB3+2 . »CC3*UO/VO+DD3*WO/VO 

GlWW=EEl+DDl*U0/V0+2 . »FF1«W0/V0 

G2WW=EE2+DD2*U0/V0+2 . *FF2*W0/V0 

G3WW=EE3+DD3*U0/V0+2 . *FF3*W0/V0 

D(1)=-GG1 

D(2)=-GG2 

D(3)=-GG3 

C(l,l)=GlWV 

C(1,2)=G1WU 

C(1,3)=G1WW 

C(2,l)=G2WV 

C(2,2)=G2WU 

C(2,3)=G2WW 

C(3,1)=G3WV 

C(3,2)=G3WU 

C(3,3)=G3WW 

DO 205 11=1,3 

DO 205 JJ=1,3 

CC(II,JJ)=C(II,JJ) 

205 CONTINUE 

CALL DLEQD(CC,D,3,1,3,3,DET) 
DELV=D ( 1 ) 

DELU=D(2) 

DELW=D(3) 

WV=VO+DELV 

WU=UO+DELU 

WW=WO+DELW 

C 

C APPLY CONVERGENCE CRITERIA 
C 

DELV=ABS ( VO-WV ) 

DELU=ABS ( UO-WU ) 

DELW=ABS ( WO-WW ) 

EPS V=EPS*ABS ( VO ) 

EPSU=EPS*ABS(UO) 

EPSW=EPS*ABS ( WO ) 
IF(DELV.GE.EPSV)GO TO 210 
IF ( DELU. GE . EPSU ) GO TO 210 
IF ( DELW. GE . EPSW ) GO TO 210 
C 

C SOLUTION CONVERGED 
C 

WRITE(6,8)*)N 

WR=WU/UTIP 

WT=WV/UTIP 

WZ=WW/UTIP 

W=SQRT ( WT*WT+UR*UR+WZ*WZ ) 

WRITE ( 6 , 61 ) WT ,WR , WZ ,W 
GO TO 220 
C 

210 WWT=WV 
WWR=WU 

ww z=ww 
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IF(N.GE.2)G0 TO 206 

VO=WV 

UO-WU 

WO=»WW 

WWT1=WV 

WWR1=WU 

WWZ1=WW 

GO TO 200 

206 V0=(WVTCWWTl)/2, 

WO»(WWZ+V7WZl)/2. 

U0=(WWR+WWRl)/2. 

wwti=wwt 

WWR1=WWR 
WWZ1=WWZ 
200 CONTINUE 
C 

C SOLUTION FAILED TO CONVERGE - USE INITIAL VALUES OF VELOCITY 
C CALCULATED FROM LINEAR EQUATIONS FOR TURBULENCE DERIVATIONS 
C 

write(6,86)nreps 

WV=VXNIT 

WU=UINIT 

UW=WINIT 

WR=WU/UTIP 

WT=WV/UTIP 

WZ=WW/UTIP 

W=SQRT ( WT*WT+WR*WR+WZ*WZ ) 

WRITE ( 6 , 61 ) WT ,WR ,WZ ,W 
C 

C USE LINEARIZED THEORY TO CALCULATE FLUCTUATING VELOCITIES 
C 

220 FACTlM . *SQRT ( VCHAN1 ( IPROBE ,I,J) ) *SQRT( VEL1 ( I , J ) ) /SL0PE1 ( IPROBE ) 
FACT2=U . *SQRT ( VCHAN2 ( IPROBE ,I,J) ) *SQRT( VEL2 ( I , J ) ) /SL0PE2 ( IPROBE ) 
FACT3=!+ . *SQRT ( VCIIAN3 ( IPROBE , I , J ) ) *SQRT ( VEL3 ( I , J ) ) /SLOPE 3 ( IPROBE ) 
FV1=SQRT ( FV1SQ ) / ATTEN 
FV2=SQRT ( FV2SQ ) /ATTEN 
FV3=SQRT ( FV3SQ ) /ATTEN 
U1=FACT1*FV1 
U2=FACT2*FV2 
U3=FACT3*FV3 
FVEL1 ( I , J ) =U1*U1 
FVEL2 ( I , J ) =U2*U2 
FVEL3(I,J)=U3*U3 
FV1P2=FV1P2 / ( ATTEN* ATTEN ) 

FV1M2=FV1M2/ ( ATTEN*ATTEN ) 

FV1P3=FV1P3/ ( ATTEN* ATTEN ) 

FV1M3=FV1M3/ (ATTEN*ATTEN) 

FV2P3=FV2P3/ (ATTEN* ATTEN ) 

FV2M3=FV2M3/ ( ATTEN* ATTEN ) 

FV12= ( FV1P2-FV1M2 )/k. 

FV13® ( FV1P3-FV1M3 ) /U . 

FV23= ( FV2P3-FV2M3 ) A . 

FVEL12=FACT1*FACT2*FV12 

FVEL13=FACT1*FACT3*FV13 
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IF(N.GE.2)G0 TO 206 

V0=WV 

UO=WU 

W0=WW 

WWT1=WV 

WWR1=WU 

WWZ1=WW 

GO TO 200 

206 V0=(WWT+WWTl)/2. 

W0=(WWZ+WZl)/2. 

U0=(WWR+WWRl)/2. 

UWT1=WWT 
WWR1=WWR 
WWZ1=WWZ 
200 CONTINUE 
C 

C SOLUTION FAILED TO CONVERGE - USE INITIAL VALUES OF VELOCITY 
C CALCULATED FROM LINEAR EQUATIONS FOR TURBULENCE DERIVATIONS 
C 

WRITE(6,86)NREPS 

WV=VTNIT 

WU=UINIT 

WW=WINIT 

WR=WU/UTIP 

WT=WV/UTIP 

WZ=WW/UTIP 

w=sqrt ( wt»wt+ur*wr+wz*wz ) 

WRITE ( 6 , 61 )WT,WR ,WZ ,W 

c 

C USE LINEARIZED THEORY TO CALCULATE FLUCTUATING VELOCITIES 
C 

220 FACTlM . *SQRT ( VCHAN1 ( IPROBE , I , J ) ) »SQRT( VEL1 (l,J) ) /SL0PE1 ( IPROBE ) 
FACT2*1+ . *SQRT ( VCHAN2 ( IPROBE , I , J ) )»SQRT(VEL2(l , J) ) /SL0PE2 ( IPROBE ) 
FACT3=^ . *SQRT( VCHAN3( IPROBE, I , J) )*SQRT( VEL3(l , J) ) /SL0PE3( IPROBE) 
FV1=SQRT( FV1SQ) /ATTEN 
FV2=SQRT ( FV2SQ ) /ATTEN 
FV3=SQRT( FV3SQ) /ATTEN 
U1=FACT1 # FV1 
U2=FACT2*FV2 
U3=FACT3*FV3 
FVEL1 ( I , J ) ®U1*U1 
FVEL2 ( I , J ) =U2*U2 
FVEL3(I,J)=U3*U3 
FV1P2=FV1P2/ ( ATTEN*ATTEN ) 

FV1M2=FV1M2 / (ATTEN*ATTEN ) 

FV1P3=FV1P3/ ( ATTEN»ATTEN ) 

FV1M3=FV1M3/ ( ATTEN* ATTEN ) 

FV2P3=FV2P3/ (ATTEN* ATTEN ) 

FV2M3=FV2M3/ ( ATTEN*ATTEN ) 

FV12* ( FV1P2-FV1M2 )/U . 

FV13= ( FV1P3-FV1M3 ) A . 

FV23= ( FV2P3-FV2M3 ) A . 

FVEL12=FACT1*FACT2»FV12 

FVEL13=FACT1*FACT3*FV13 



FVEL23=FACT2*FACT3*FV23 
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C 

C CALCULATE FLUCTUATING VELOCITY COEFFICIENTS 
C 

AA1=SQRT(A4) 

AA2=D4/(2.*AA1) 

AA3=E4/(2.*AA1) 

aa4=aai* ( b4 / ( 2 . *A4 ) -D4»d4 / ( 8 . *A4 *a4 ) ) 
AA5-AA1* (c4/(2.*A4)-E4*e4/(8. «a4*A4 ) ) 
AA6=AA1* ( F4 / ( 2 . *A4 ) -D4*e4/( 4 . *A4*A4 ) ) 
AA7=AAI*(b4/a4-d4*i>4/ ( 4 . »a4*A4 ) ) 

AA8=AA1* ( c4/a4-e4*e4/( 4 . »A4*A4 ) ) 

AA9=AA6 

BB1=SQRT(A5) 

BB2=D5/(2.*BBI) 

BB3=E5/(2.«BB1) 

BB4=BB1* ( B5/( 2 . »A5 ) -D5*D5/ ( 8 . *A5*A5 ) ) 
BB5=BB1*(C5/(2.*A5)-E5*E5/(8.*A5*A5)) 
BB6=BB1* ( F5 / ( 2 . »A5 ) -D5*E5/ ( 4 . «A5*A5 ) ) 
BB7=BB1* ( B5 /A5-B5«D5 / ( 4 . *A5*A5 ) ) 

BB8=BB1* ( C5/A5-E5*E5/( 4 . *A5*A5 ) ) 

BB9=BB6 

CC1=SQRT(A6) 

CC2=D6/(2.*CCl) 

CC3=E6/(2.*CC1) 

CC4=CC1« ( B6/ ( 2 . *A6 ) -D6*d6 / ( 8 . *a6«a 6 ) ) 
CC5=CC1* ( C6/ ( 2 . *A6 ) -E6*E6/ ( 8 . «A6«A 6 ) ) 
cc6=cci* ( f6 / ( 2 . *A6 ) -d6*e6 / ( 4 . »A 6 *A 6 ) ) 
CC7=CC1* ( B6/A6-D6*D6/ ( 4 . *a6«A6 ) ) 

CC8=CC1 # ( C6 /A6-e6*e6 / ( 4 . »a6*A6 ) ) 

CC9=CC6 

AAA1=AA1 

AAA2=BB1 

AAA3=CC1 

BBBl=AA2+AA7* T -m /UT+AA6*WZ /WT 

BBB2=BB2+BB7*WR/WT+BB6*WZ/WT 

BBB3=CC2+CC7*WR/WT+CC6*WZ/WT 

CCC1=AA3+AA8*UZ/WT+AA6*WR /WT 

CCC2=BB3+BB8*WZ/WT+BB6*WR/WT 

CCC3=CC3+CC8*VZ/WT+CC6*WR/WT 

A(l,l)=AAAJL*AAAl 

A ( 1 , 2 ) =BBB1*BBB1 

A( 1 , 3 ) =CCC1*CCC1 

A( 1 ,4 )=2 . *AAA1*BBB1 

A ( 1 , 5 ) =2 . *AAA1*CCC1 

A ( 1 , 6 ) =2 . *BBB1* CCC1 

A(2,1)=AAA2«AAA2 

A(2,2)=BBB2«BBB2 

A( 2 , 3 ) =CCC2*CCC2 

A( 2 , 4 )=2 . *AAA2«BBB2 

A(2,5)=2.*AAA2*CCC2 

A(2,6)=2.*BBB2*CCC2 

A(3,1)=AAA3*AAA3 

A ( 3 , 2 ) =BBB 3*BBB3 
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A(3,3)=CCC3»CCC3 

A(3,U)=2.«AAA3*BBB3 

A(3,5)= s 2.*AAA3»CCC3 

A(3,6)=2.*BBB3»CCC3 

A(U,l)=AAAl*AAA2 

A(U,2)=BBB1«BBB2 

A(4,3)=CCC1»CCC2 

A(4,4) =AAA1*BBB2+AAA2»BBB1 

A( 4 , 5 ) =AAA1*CCC2+AAA2*CCC1 

A(4, 6)=BBB1*CCC2+CCC1*BBB2 

A(5,1)=AAA1»AAA3 

A ( 5 » 2 ) =BBB1*BBB3 

A ( 5 , 3 ) =CCC1*CCC3 

A( 5 ,4)=AAA1»BBB3+AAA3*BBB1 

A( 5 , 5 )=AAA1*CCC3+AAA3»CCC1 

A(5»6) =BBB1 # CCC3+CCC1*BBB3 

A(6,l)=AAA2*AAA3 

A(6,2)=BBB2*BBB3 

A(6,3)=CCC2*CCC3 

A( 6 ,4 )=AAA2«BBB3+AAA3*BBB2 

A ( 6 , 5 )=AAA2*CCC3+AAA3*CCC2 

A( 6 ,6 )=BBB2«CCC3+CCC2«BBB3 

b(i)=fveli(i,j) 

B(2)=FVEL2(I,J) 

B(3)=FVEL3(I,J) 

B(4)=FVEL12 
B(5)=FVEL13 
B(6)=FVEL23 
DO 113 11=1,6 
DO 113 JJ=1,6 
AA(II,JJ)=A(II,JJ) 

113 CONTINUE 

CALL DLEQD(AA,B,6,1,6,6,DET) 

B(l)=LSqRT(B(l)) 

B(2)=DSQRT(B(2)) 

B(3)=DSQRT(B(3)) 

C 

C CORRECT FLUCTUATING VELOCITIES FOR HIGH TURBULENCE INTENSITIES 
C 

ARQ1=1 . +FACTOR* ( B ( 1 ) / ( W»UTIP ) ) **2 . 

ARG2=1 . +FACT0R# ( B ( 2 ) / ( W*UTIP ) ) *»2 . 

ARG3=1 . +FACT0R* ( B ( 3 ) / ( W«UTIP ))**2. 

B(1)=B(1)/SQRT(ARG1) 

B(2)=B(2)/SQRT(ARG2) 

B(3)=B(3)/SQRT(ARG3) 

C 

FWT=B(1)/(W«UTIP) 

fwr=b(2)/(w*utip) 

FWZ=B ( 3 ) / ( w*utip ) 

WRITE (6,62) FWT ,FWR ,FWZ 
B ( 4 ) =B ( k ) / ( W*W*UTIP*UTIP ) 

B ( 5 ) =B ( 5 ) / ( W*W*UTIP*UTIP ) 

B ( 6 ) =B ( 6 ) / ( W*W»UTIP*UTIP ) 

write(6,63)b(4),b(5),b(6) 



QSOJ^^'ft^FWR+FfrZ^FWZ+FWT^FWT 

WRITE(6,67)QSqD 

RATI01=B(lt)/QSQD 

RATI02=B(5)/QSQD 

RATI03=B(6)/QSQD 

WRITE ( 6 , 68 )RATI01 ,RATI02 ,RATI03 

WRITE(6,99) 

110 CONTINUE 
115 CONTINUE 

FORMAT STATEMENTS 

5 FORMAT (9X, II) 

6 FORMAT (3F10. 5) 

10 F0RMAT(8X,I2,F10.5) 

15 FORMAT ( ItFlO . 5 ) 

20 F0RMAT(9X,I1,8X,I2) 

25 FORMAT (10F7. 3) 

32 FORMAT (6F10. 5) 
kO FORMAT (9X, II, 9X, II) 

k5 FORMAT (3F10. 6 ) 

U 6 FORMAT (U FI 0 . 6 ) 
hj F0RMAT(UF10.6,9X,I1) 

6 1 FORMAT (1 OX, 'MFAN VELOCITIES :',UE15.5) 

62 FORMAT (1 OX,' TURBULENCE INTENSITIES :',3E15.5) 

63 FORMAT (1 OX,' TURBULENCE STRESSES ;',3E15.5) 

67 FORMAT (1 OX, 'TOTAL TURBULENCE ENERGY :',E 15 . 5 ) 

68 FORMAT (1 OX, ' STRESS/ENERGY RATIOS :',3E15.5) 

70 F0RJ1AT(3F10.5) 

75 FORMAT(6F10.6) 

82 FORMAT(5X,I5,F10.5) 

8 U FORMAT ( 5X , ' MEAN VELOCITIES CO T TVFRGET) IN' ,15, 'ITERATIONS' ) 

86 FORMAT (5X,' MEAN VELOCITIES DID NOT CONVERGE AFTER ', 15 ,' ITERATIONS ' 

1 /) 

99 FORMAT ( IX , / / ) 

STOP 
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